{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## *This notebook reproduces entanglement spectrum results in Li and Haldane, PRL 101, 010504 (2008)*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Load QHED code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "qhed_dir = \"../\"\n",
    "# CG_dir = joinpath(qhed_dir, \"CG\") # specify directory for precomputed Clebsch-Gordan coefficients\n",
    "\n",
    "include(joinpath(qhed_dir, \"src\", \"QHED.jl\"));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set matplotlib defaults"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "using PyPlot\n",
    "\n",
    "PyPlot.matplotlib[:rc](\"font\", size=14)\n",
    "PyPlot.matplotlib[:rc](\"text\", usetex=true)\n",
    "PyPlot.matplotlib[:rc](\"xtick\", direction=\"in\", bottom=\"on\", top=\"on\")\n",
    "PyPlot.matplotlib[:rc](\"xtick.minor\", visible=\"on\")\n",
    "PyPlot.matplotlib[:rc](\"ytick\", direction=\"in\", left=\"on\", right=\"on\")\n",
    "PyPlot.matplotlib[:rc](\"ytick.minor\", visible=\"on\")\n",
    "\n",
    "# bblue = \"#0C2340\"\n",
    "# borange = \"#FC4C02\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Set parameters and build setup\n",
    "$N$ = number of electrons\n",
    "\n",
    "$N_\\phi$ = number of flux quanta penetrating sphere\n",
    "\n",
    "For a LLL problem at filling factor $\\nu$, we will in general have $N_\\phi = \\nu^{-1}N - S$, where $S$ is the \"shift\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Building states list for N = 16, Nphi = 29, Lz = 0 ...\n",
      "Finished building states list of size 2374753\n",
      "elapsed time: 355.810371755 seconds\n",
      "\n",
      "Building L^2 matrix of size 2374753 x 2374753 ...\n",
      "Working on LpLm, column j = 1: 1/2374753 = 4.2109642560721055e-5%\n",
      "Working on LpLm, column j = 118738: 118738/2374753 = 5.0000147383748965%\n",
      "Working on LpLm, column j = 237476: 237476/2374753 = 10.000029476749793%\n",
      "Working on LpLm, column j = 356214: 356214/2374753 = 15.000044215124689%\n",
      "Working on LpLm, column j = 474952: 474952/2374753 = 20.000058953499586%\n",
      "Working on LpLm, column j = 593690: 593690/2374753 = 25.00007369187448%\n",
      "Working on LpLm, column j = 712428: 712428/2374753 = 30.000088430249377%\n",
      "Working on LpLm, column j = 831166: 831166/2374753 = 35.00010316862427%\n",
      "Working on LpLm, column j = 949904: 949904/2374753 = 40.00011790699917%\n",
      "Working on LpLm, column j = 1068642: 1068642/2374753 = 45.000132645374066%\n",
      "Working on LpLm, column j = 1187380: 1187380/2374753 = 50.00014738374896%\n",
      "Working on LpLm, column j = 1306118: 1306118/2374753 = 55.00016212212386%\n",
      "Working on LpLm, column j = 1424856: 1424856/2374753 = 60.000176860498755%\n",
      "Working on LpLm, column j = 1543594: 1543594/2374753 = 65.00019159887366%\n",
      "Working on LpLm, column j = 1662332: 1662332/2374753 = 70.00020633724854%\n",
      "Working on LpLm, column j = 1781070: 1781070/2374753 = 75.00022107562344%\n",
      "Working on LpLm, column j = 1899808: 1899808/2374753 = 80.00023581399834%\n",
      "Working on LpLm, column j = 2018546: 2018546/2374753 = 85.00025055237323%\n",
      "Working on LpLm, column j = 2137284: 2137284/2374753 = 90.00026529074813%\n",
      "Working on LpLm, column j = 2256022: 2256022/2374753 = 95.00028002912303%\n",
      "Working on Lz^2, column j = 1: 1/2374753 = 4.2109642560721055e-5%\n",
      "Working on Lz^2, column j = 237476: 237476/2374753 = 10.000029476749793%\n",
      "Working on Lz^2, column j = 474952: 474952/2374753 = 20.000058953499586%\n",
      "Working on Lz^2, column j = 712428: 712428/2374753 = 30.000088430249377%\n",
      "Working on Lz^2, column j = 949904: 949904/2374753 = 40.00011790699917%\n",
      "Working on Lz^2, column j = 1187380: 1187380/2374753 = 50.00014738374896%\n",
      "Working on Lz^2, column j = 1424856: 1424856/2374753 = 60.000176860498755%\n",
      "Working on Lz^2, column j = 1662332: 1662332/2374753 = 70.00020633724854%\n",
      "Working on Lz^2, column j = 1899808: 1899808/2374753 = 80.00023581399834%\n",
      "Working on Lz^2, column j = 2137284: 2137284/2374753 = 90.00026529074813%\n",
      "elapsed time: 464.154080105 seconds\n",
      "\n"
     ]
    }
   ],
   "source": [
    "N = 16\n",
    "Nphi = 2N - 3 # choose the Moore-Read shift\n",
    "\n",
    "setup = HaldaneSphereSetupLLL(N, Nphi, 0);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate Coulomb pseudopotentials for the first excited ($n = 1$) Landau level"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "SphericalPseudopotentials(14.5, 13.5, [29.0, 28.0, 27.0, 26.0, 25.0, 24.0, 23.0, 22.0, 21.0, 20.0  …  9.0, 8.0, 7.0, 6.0, 5.0, 4.0, 3.0, 2.0, 1.0, 0.0], [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0  …  20.0, 21.0, 22.0, 23.0, 24.0, 25.0, 26.0, 27.0, 28.0, 29.0], [0.636754, 0.43937, 0.486364, 0.342107, 0.288476, 0.256355, 0.234165, 0.217661, 0.204809, 0.194478  …  0.146247, 0.144548, 0.143087, 0.141844, 0.140804, 0.139955, 0.139287, 0.138792, 0.138465, 0.138302])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "pp = spherical_Coulomb_pseudopotentials(get_ell(setup), 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Build Hamiltonian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Building CG table for ell1 = 14.5, ell2 = 14.5:\n",
      "Working on alpha1 = 1: 1/30\n",
      "Working on alpha1 = 2: 2/30\n",
      "Working on alpha1 = 3: 3/30\n",
      "Working on alpha1 = 4: 4/30\n",
      "Working on alpha1 = 5: 5/30\n",
      "Working on alpha1 = 6: 6/30\n",
      "Working on alpha1 = 7: 7/30\n",
      "Working on alpha1 = 8: 8/30\n",
      "Working on alpha1 = 9: 9/30\n",
      "Working on alpha1 = 10: 10/30\n",
      "Working on alpha1 = 11: 11/30\n",
      "Working on alpha1 = 12: 12/30\n",
      "Working on alpha1 = 13: 13/30\n",
      "Working on alpha1 = 14: 14/30\n",
      "Working on alpha1 = 15: 15/30\n",
      "Working on alpha1 = 16: 16/30\n",
      "Working on alpha1 = 17: 17/30\n",
      "Working on alpha1 = 18: 18/30\n",
      "Working on alpha1 = 19: 19/30\n",
      "Working on alpha1 = 20: 20/30\n",
      "Working on alpha1 = 21: 21/30\n",
      "Working on alpha1 = 22: 22/30\n",
      "Working on alpha1 = 23: 23/30\n",
      "Working on alpha1 = 24: 24/30\n",
      "Working on alpha1 = 25: 25/30\n",
      "Working on alpha1 = 26: 26/30\n",
      "Working on alpha1 = 27: 27/30\n",
      "Working on alpha1 = 28: 28/30\n",
      "Working on alpha1 = 29: 29/30\n",
      "Working on alpha1 = 30: 30/30\n",
      "elapsed time: 315.37450272 seconds\n",
      "\n",
      "Building 2-body interaction matrix ...\n",
      "Working on alpha1 = 1: 1/30\n",
      "Working on alpha1 = 2: 2/30\n",
      "Working on alpha1 = 3: 3/30\n",
      "Working on alpha1 = 4: 4/30\n",
      "Working on alpha1 = 5: 5/30\n",
      "Working on alpha1 = 6: 6/30\n",
      "Working on alpha1 = 7: 7/30\n",
      "Working on alpha1 = 8: 8/30\n",
      "Working on alpha1 = 9: 9/30\n",
      "Working on alpha1 = 10: 10/30\n",
      "Working on alpha1 = 11: 11/30\n",
      "Working on alpha1 = 12: 12/30\n",
      "Working on alpha1 = 13: 13/30\n",
      "Working on alpha1 = 14: 14/30\n",
      "Working on alpha1 = 15: 15/30\n",
      "Working on alpha1 = 16: 16/30\n",
      "Working on alpha1 = 17: 17/30\n",
      "Working on alpha1 = 18: 18/30\n",
      "Working on alpha1 = 19: 19/30\n",
      "Working on alpha1 = 20: 20/30\n",
      "Working on alpha1 = 21: 21/30\n",
      "Working on alpha1 = 22: 22/30\n",
      "Working on alpha1 = 23: 23/30\n",
      "Working on alpha1 = 24: 24/30\n",
      "Working on alpha1 = 25: 25/30\n",
      "Working on alpha1 = 26: 26/30\n",
      "Working on alpha1 = 27: 27/30\n",
      "Working on alpha1 = 28: 28/30\n",
      "Working on alpha1 = 29: 29/30\n",
      "Working on alpha1 = 30: 30/30\n",
      "elapsed time: 33.806678304 seconds\n",
      "\n",
      "Building Hami matrix of size 2374753 x 2374753 ...\n",
      "Working on H, column j = 1: 1/2374753 = 4.2109642560721055e-5%\n",
      "Working on H, column j = 118738: 118738/2374753 = 5.0000147383748965%\n",
      "Working on H, column j = 237476: 237476/2374753 = 10.000029476749793%\n",
      "Working on H, column j = 356214: 356214/2374753 = 15.000044215124689%\n",
      "Working on H, column j = 474952: 474952/2374753 = 20.000058953499586%\n",
      "Working on H, column j = 593690: 593690/2374753 = 25.00007369187448%\n",
      "Working on H, column j = 712428: 712428/2374753 = 30.000088430249377%\n",
      "Working on H, column j = 831166: 831166/2374753 = 35.00010316862427%\n",
      "Working on H, column j = 949904: 949904/2374753 = 40.00011790699917%\n",
      "Working on H, column j = 1068642: 1068642/2374753 = 45.000132645374066%\n",
      "Working on H, column j = 1187380: 1187380/2374753 = 50.00014738374896%\n",
      "Working on H, column j = 1306118: 1306118/2374753 = 55.00016212212386%\n",
      "Working on H, column j = 1424856: 1424856/2374753 = 60.000176860498755%\n",
      "Working on H, column j = 1543594: 1543594/2374753 = 65.00019159887366%\n",
      "Working on H, column j = 1662332: 1662332/2374753 = 70.00020633724854%\n",
      "Working on H, column j = 1781070: 1781070/2374753 = 75.00022107562344%\n",
      "Working on H, column j = 1899808: 1899808/2374753 = 80.00023581399834%\n",
      "Working on H, column j = 2018546: 2018546/2374753 = 85.00025055237323%\n",
      "Working on H, column j = 2137284: 2137284/2374753 = 90.00026529074813%\n",
      "Working on H, column j = 2256022: 2256022/2374753 = 95.00028002912303%\n",
      "elapsed time: 2300.086236136 seconds\n",
      "\n",
      "2649.575526 seconds (10.88 G allocations: 1.949 TiB, 14.13% gc time)\n"
     ]
    }
   ],
   "source": [
    "@time Hami = HaldaneSphereHamiLLL(setup, pp, build_CG_table(get_ell(setup)));\n",
    "# @time Hami = HaldaneSphereHamiLLL(setup, pp, load_CG_table(get_ell(setup), CG_dir));"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Diagonalize Hamiltonian to find ground state wave function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Diagonalizing Hami of size 2374753 x 2374753 with eigs, nev = 1 ...\n",
      "elapsed time: 531.746492544 seconds\n",
      "\n"
     ]
    }
   ],
   "source": [
    "eigs!(Hami, 1)\n",
    "psi = Hami.eigenstates[:, 1];"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Calculate and plot entanglement spectra for different bipartition configurations (cf. Fig. 2 of Li and Haldane)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculating entanglement spectrum for NA = 8, LzA = [49.0, 50.0, 51.0, 52.0, 53.0, 54.0, 55.0, 56.0, 57.0, 58.0, 59.0, 60.0, 61.0, 62.0, 63.0, 64.0, 65.0, 66.0, 67.0] ...\n",
      "elapsed time: 696.692456998 seconds\n",
      "\n"
     ]
    }
   ],
   "source": [
    "subsystemA = collect(16:get_Norb(setup))\n",
    "NA = 8\n",
    "LzAvec = collect(49.:67.)\n",
    "\n",
    "ent_spec = entanglement_spectrum(psi, setup.states, subsystemA, NA, LzAvec);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjIAAAHUCAYAAAAgOcJbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3c+LI2ma4PnH0yvbi5kqD4VHsVBFBRNpPhR0sCyNPAR76aa3Q151GGZhu6XMvyClXXxYltxGSodmM3O3KV85e1gGRCPFfSBCXnPLw4wU0HQflu0IiT4Mw1CMrAuimL5UuJsri9nyrvGyPUSayuUuk0zSa7LnffX9QBC4XDJ/7bXX9D72/twKwzAUAAAAC72XdQIAAACWRSADAACsRSADAACsRSADAACsRSADIHPValV6vV7WycAcXCdoRCADIFO9Xk/a7bYEQZB1UjAD1yk9vV5P7t+/L/fv35fDw8O1//3T09Px369Wq2v/+6sikAGQqW63KyIivu9nnBLMwnVK14cffigXFxfjfF6nWq0mFxcX0ul01v63TfhG1gmAfQaDgZycnMwt9NH7Ip7nSaPRSDt5zmi329LtdscVx0cffSS1Wm3qeweDgbRaLREROT8/l0KhEPteTU5PT6VarcrZ2ZkMh8Osk7OSefdFtVqV/f19KZVK4nnexGeOj48ln8+vM7kL0X6der3eRL7ncjk5Pj6WXC6XYarMiM4tl8tJEASyv78/9d72fV8ajcb4e2CZY1grxFqUSqUwn8+HIhKKSHhxcRH7vlwuF4pI6HleWKlU1pzSeBcXF2Gr1QpFJCwWizPf22q1wnw+Hw6Hw/FrnU4nrNVqRtLiQn7OUqvVwk6nM/55OByGuVwu9DzvzntbrVboed5EHpRKpbBUKhlLTxr5PRwOw0ajEYZhGBaLRaPpXaek90WpVBrn381/UR6YsInXqdVqha1Wa+K1brcb5vP5taajVCqFnueN8z7Kq1KpFBaLxTCfz4fFYjHs9/t3Ptvtdqdeg06nc+c8Go3G1GvQ7/enlr9FjhGXDu0IZNYoKmgiMrNC73Q66r4sopux1WqFuVxu5hd2v98Pc7ncnS9Rz/OMfrnYnJ+z9Pv9qefT7XZDEZn4orm4uAhFZCLoufl6t9s1mi6T+X3zPKIK2DaL3Bc3g4zovdMqtVVt0nW6uLiIzfNGo2E0SEyi3+/fuUdvigLe2/drXACRy+Wm3sO5XO7OMeICmUWOQSCDuRqNxvjJOpfLxb6v1Wql8gVnyrwv7LgnPNMtIq7k522VSiX2STp6ir753rgnb8/z5racLcJkfrdarYkv11qtNvOYNkgSyKzDJl2nWRVvXMWeprhA5abb93AYTj+PTqdz532RaS2u08530WPYGsgw2HeNhsOheJ4nH374oQRBEDuNcTgcqu4vn6XX64nv+1Iul+/8bjgcTu2/XZar+fn69Wv54IMPps4OicZWRONmer2e5HK5qWMB8vm80amypvI7CAJptVrS7XalXq9LvV6XwWDAbBhDNuk6zTq/8/PztY+RiQbqFovF2Pfs7e2JyPxB091uNzb9nuclurdNHMMGBDIZqNfrIiKxA181fVEsKhpwN+tGNs21/Nzb25MgCOT8/Hzue2d9GSb9wlzUqvldr9fl5cuX0mg0xv+iKZ+2XSvNNuE65fN58X1fDg4O7qSp1WqtfSpxr9cTz/NmBlDR/Rjdn/OONc3+/r4EQSCDwSD1Y9iAQGZNBoOBHBwciMi7SDh6Wr5980U3pa1uRvinp6dyeno6fpozyeX87HQ60u/3p34BRV+CcV9ON0Vfpje/qKKZM4syld/tdlsODw/vfNHfbmlyVRAE4/tB831hy3XyPE9KpZIMBgP54IMPxt8/7XZbPM9b6IFq2Xsj4vu+BEEw82+enZ2JiEipVJrbWpTkQWbee0wcwwYEMmvS6/UmCvjx8bGIyMT05Gnvs43v+5LL5aTdbkulUpFarSaNRkP29/fHTwAmuJyfuVxuapN/VCHdnDY5q2tgWmUTdf0tWhGtmt9Rd2O1Wr2zTsbZ2dm4Qq/X6840d09zcnIy0cIhIkYXQNvE69TpdKRUKkkQBHJ4eCgHBweyt7e38FIPy94bNz8vIjMDxHq9Lvl8PtF6LUEQxLbaRK/P+z41cQwrZD1IZ1NMG0AlX89euGmZ6cnRALxl/i0zq2XWoMbonG5PhwzDd4NPTU2/TjM/51l3fkfy+fydada1Wm3mYF8RmbgWw+FwqYHPWeb3sjTdF2EYxub77Wu0Cm3XaV3X4OLiYjyIVr6eUn5z6Ycklr03ItH0+ml/t9/vjyc7TLtXpw2ylRnT+aNBvDfLzbTBvosew9bBviyIl6FKpSLtdnvi6WiZ6Pjm050GQRDIhx9+eOf1fD4vp6enqS1UtWx+np2dSalUSvx3ssjv09NTGQwGMhwOJ/Lu+PhYTk9PpdfrTZyD7/vieZ74vj/xRJakSyopU+U3Ldrui7jWs3w+L41GQyqVSip/N8vrtI5r0Ov1pFqtSqfTkcFgIPV6fdx1FtdFO82q90bUIjNtQsODBw+k2+0avf/wW3QtrUFcf/TtwXi+76/UR6vJtEDFVP+66fzU3kfs+76cnJxM/SLM5XLS7Xbl5ORkXDkFQSBnZ2fj964aNG5i+V2nvb09I2NONvE6RV1h3W5XisWi1Gq1cfASBMHU2ZNppSMaH3Oz6zD6V6vVjAYx0XfWvAHDaR9DCwKZNYjrj749GG8wGFg3nuO2uKnAIu+eSkTeTS9exSblp4hIuVyWTqcTey7FYlE6nY602205PT2VdrsttVpt/EX15MmTlf5+FvkdBIGq1h0T4s5n2qDsZWzafSHyboDu8fHxRJCQz+dlOBxKsViUwWCwllk5UWuMycApl8vNfcia95Bi4hg2IJBZg1nNm9H0wHa7La9evbJqvZNpnjx5EvuF/fbtWxFZvQnXZH5GXTBalctlOT4+nlvxeJ4ntVpt/E/kt+e26hfVustvr9eTp0+fytOnT1c+lhblclnu378/9d6IXtN0X9jidpfqTdEaKuuYYZVk/ZhFRa1K00T7Xc0rMyaOYQMCmTWYVZFE/eI3uwYWVa/Xx1uwL/rP9MyDaAZGml/YJvKz3W5Lu92WRqMhvu+P+9WTWFd+1+t1OTw8vPNF3W63E6U1Gi+wqrTL721R8/yqNN0Xs4LKaOG2VQPOdV+nJNZxDWZ1jTx58mQtLQ7RwpQmg4KPPvoo9j4PgiDR3zNxDCtkPdrYdcPhcO6MhGiZ+XXvC7KsWbMzhsNh7B4/+Xx+6j4tccvxxx1/1fwsFovj9EXHuri4ULWHTKfTiU3/zWXFu91umMvl7sy2iPb+uW3RmRlZld8sNv1b1az7olarTS3n0Z5Y02YRrfu+sNHNe3ma27P8Zrm4uFhq1lL0nbfKFhTTZgtFx52WpmnbvUybtbToMWydtUSLTMoajcbciDdq9rWl33rWqrOe50mlUrnzRB31VT979mzi9YODA7l//37i1pBV87Narcrh4eGd3+W+3t5ewyqXvu/Lxx9/LM+fP5eDg4OJf/v7+xNpnJZv0aJrt/NaRMbrbCR9Knex/KZl1n1xfHw8tXWsXq9LLpe7c7+s+74QeVeWTk9P5ezsbPxPu5urDd/WbrfH+ZtEdI8t2hUV5ZPJ9YBEZLyg3/PnzydeHwwG4vt+olZLE8ewQtaRlKsqlcp4HY9cLjc3Wl/35maLqtVqYbFYvLNWQ7FYnLpBWrFYDCuVStjv98NWqxV6njf1yalUKoWe5yV6mjSRnzeLfL/fn3hSkZj1WNYt2iE57t/tcyuVSmGn05nI67hN6yqVSqI1Okzl98XFRVir1cJOpxN2Op2Jlomo1SUqQ5VKZfw0GP0u2rSw1WqpWqMmssh90e12w1KpFHa73fGTbz6fn1rm1n1fDIfDid8Nh0OV+T1Np9MZl5V+vz/O20VbnqK8TKpWq03cq57nhaVSaeH1a8IwviUkaimOykG/3w/z+fzU+zduk8xFjmFriwyBDFITVaydTmdugGBqQbBZoh2Bp/3NVqulPpicJeqKarVac/M6+rJfh3w+PxEsDofDicqi0+mMFy/r9/vjyqfb7d4JLDudjtXXKAzfVSrRtUrSjbGO+yIMf7tYZavVGt+ztul2u2Gj0Uj0faPNvAAiybnN2+07yTFsDWRYEA+pyefziWZHDAaDtQw48zxvvCHjzebmwWAgnU4n0bLhWi2yoN/5+fla8vvs7EyCIJgoA9GA12gLi2h6aJSeae+NlEolKZfLMhgMrJ11k8vlEl+rdd0XQRA40c1QLBad7d40cW4u5w9jZJC558+fr+0GixaPOzs7k36/L+12W16/fj1zu3vXrKuCfPXq1dQ83dvbG0/9FFl8FpuWDQvTtq77YtasFsAGBDLI1NnZWexgvTR4njceGFmv16VSqaS2NLxG66ycCoXC1Epy1ZVlnZguOsc674t8Pj8e7B7xfV/NppDAPAQyyFSpVMqkYnr9+vVGVIi3vXjxYmL37DSVSiXZ29ubOssqSfAYLfseabfbUiwWre1WWsS674tOpzNuqTw7OxPf9xfqrgSyxBgZYIOsu/Wp3+/LycnJeFuK4XAo/X5fRN4tItZoNMYL91Wr1YnKu9VqjVsFzs/PZTgcjldQhVkuj5+wxYsXL+TFixfy5MmTpcr5KovbnZ6eysnJiYjI1A1/tdsKwzDMOhEAAADLUNe1NGtZ9V6vN35yOzw8lHa7vebUAQAATVR1LQ0GA3n69OnU5u9oZdhoimAQBPLBBx9Iv9+XVqu17qQCAAAFVLTI+L4v5XJZnj9/HrsBWKvVmhikGC3rnXQDPQAA4B4VgYznedLpdKTRaMSu5fHixQs5PT2deO3JkyciIkwTBABgQ6kIZJLY29uTt2/fZp0MAACgiKoxMrPcXAk0EnUpRS0zAABgs1jTIjNNq9XamAWyAADAXda0yNwWrT4ZLa4V5ze/+Y387Gc/k/fff1+2trbGr+/s7MjOzk7ayQQAAAlcXV3J1dXV+OcwDOXXv/61PHr0SN57L77dxcpAJggCqdfriTb6+9nPfrbSvi4AACA7b968ke9///uxv7cykCmXy9LtdhMtx/z++++LiMjf/M3fyHe/+93x64u2yBQKBXn16tXiiVV8HBPHGI1G8vDhQ3nz5o3s7u5mnh5Tx9GUFlN5nMU5/ed/+C+xv/uDP/gD+au/+qupv/tHvzP51fRff/Zvkifwhn/3xY/mpuWrr76SH/zgB/LTn/5Uvv3tb89Ny9tf/ir27/3wRz+Sf/tvpqf1wbe+OfGziXOaxdbrZPI4cWzNG235a/o63W6R+eqrr+Tx48dT78ubrAtkqtXqePfiyGAwiB0nE3Unffe7350Z0c2zvb29ckWt7Tim0iIisru7qyY9mvJGUx5ncU7/zadfxv/yv2/If/t//T9Tf/Wz//OfTfz83s4/Spy+m26mc1Za/sn/+q/lsPXvEqVl1rm/f30lH3zvv0qUNhPnNMsi1+kbMyrrWce5XVmbOicTx5kVgLz3O9+Ub3xz+t+4fU6myvB/aPwPsYc5ODiIHSZxMz2D/+OfT33PV6Nfyj/9p/vyH//jUL69+607v9+9FUSbuk4mzmmW0WgkIjIxLGQaqwKZdrst5XJ5ImjxfV983099wO/R0ZFzxzGVFlNczBtNebxIWmZVAh//T/8i9vdJv6AW9e//9/gn7VarLdXqejfDjGPr9X78v814sv5nP479/e3KOu46jUZfyfe+9135T//p72V3d/bTtSmmzsmUWffGv/gfK4nune/cCkgiv/Obf5Df/H8jefCt37kTtGSlUqmu7W+p2zRyf39fisXinW0Hop1yDw8PJ17vdrvSaDRiA5mf//zn42b5VVpkMN1oNJJ79+7J5eWlsZYHTMoijx/Negqd4XYlMCsgmiWNgOgXMV1C855mb1ceps7JRB5rSsssi5ZhE+dl6pzSzhsTsshfET3lRkWLTBAEcnJyMm5dabfbcn5+LoVCYbwtQblcliAIpq7iO6s1JhoHwwyldOzs7Mhnn31G/qbI5jxOq4VmGU/+/GXs7x7+z/9K/rt/+f9O/d3tL92ZT/oz3D7OX9f+cKnjpJGWWS1eJixahk2Um9d/9nTlY5iUZlBv83eECepaZEyjxQBYnKaWFBE3n9BNHMeG1oKsaGtVdLEFjhYZAGppakkRMdPyYKrVIe3WiyxoC1w10XSOplrgXKPnCgFAikxVSKaOY6JryVT3iakKUlNApC3g1NbV5RICGQDqaaqUTFXWv3/6l0sd52bwMGvcT9JjmGQqIIobmD3PzYHZmlpSRMxcK033gYie9Oi60gAwhYlKyVQAQvN++kxU+ppah1ylJa90pAIAUqYtANHyNCtipptLGxevt6lzcq2bi0AGABZgKgDRNMXYRDeXSa5VtCJ6Wi9E9HVJrkpPzgKABUxVSCa6PrRVSKaCvLgVbBdhKhgy1UVl4jiaWvFE9HTfEcgAwAJcHGvjYguIqSDP1HUycRxTAYCpgEhLGd6YQKZQKMj29rYcHR2p2g8FcJmWJzYRfVOVTTBVIWmr9DGdqftJUzfXNM1mU5rNplxfXyd6v+6zMejVq1es7AusmaaKzcVuGO0VEn5L02Bf7aIGh2hl33m4CwAgA5qCEG1dS5pa8kzRnLZlaRmz417OAlBDyxedRpoqa22tVZpaHjSVYU1p0YRABkBqXHwKNcVEZa0pGBJxs6LVNEtN2/2kJeDUlSsAkBIXK1ktFYlG2q63pmulLQBelc5UAcANmr54tY0n0cRUZW1ipWGtle4qtE391xIsunelAaxMU+AgYuaL19SXt4nF2kTMBERaKhLTTKw0rK0MuzhrSUuwqCMVAIx88Wp7YnORqTzWtDGiqYBIU2BFGd4cBDKAEppaHUzRVNG6uCCeqdlGpp6sTR3Hxe47TeXGNQQyAO7QtoS5iQpS2xRjTa0X2rphTHTfacpfbVzLGwIZQAkTXy6admY2SVNFq6kS0NbKpGmDRW1cvL+1XCe9Vx3YMCZubm1f5JpadkxV+qbyWNMmgqZo2mBRW3Cm6f7+xS9/tdTnbreUaeku03UXAHCKpopWW9eSCa5No9VIS2Vtkmv3gp5vGQCIQUVrB65TPC3dMC7amBwqFAqyvb093lUTQPo0Nctrq2Q1pUfToGxtNHWPmqJ9Vliz2ZRmsynX19eJ3r8VhmGYcpoyFW0Dfnl5Kbu7u1knB9gojz79cqnPaW3C1sbUWAdT10nTYF9tLSAu3gtp53HS+tu98BmAc7RVSlqYGuugbfaTpoGxiKclj3WkAoCTtFWQJmhqdTBF2+BPTXljiqauRNfoveoArKetgjRB0xRjV7mYN5qDLNuRswCw4VxsLXCxVQfTccUApMbFzQhdPCdtlbeLO0UjPbpKLwCnaNuM0ARNaTE1Bklb64WmPNaUN5rSoonbZwcASploMTA1Bklb64WJaeUurv+iKS2aEMgAuIMnPyzDVPDg4iBxpIdvHQB38OSXPhPdQprG2Wjj4j5UmtKiCYEM4BBaUuxhotVB23UzFTxoWkJfUx5rSosm5AqwIk0LpLn4FKoNwWL6XNxbC+nhzgJW5OICaVS68VwMFjW1gIjouhcIXPUjpwGHaKocMZuJis5UJevi4FoXHzAwHYEMsCJNC6S5+hSo6alYU7DoaiWrrYUIurn5rTdFoVCQ7e1tOTo6kqOjo6yTA4e4uOibKaYCEE0VtovXSVNwJmKmhUjTAwYW02w2pdlsyvX1daL3b4VhGKacpkyNRiO5d++eXF5eyu7ubtbJAaxgKgB59OmXSx3ndgBi6jiu0dRSJaKv3MBuSetv9x4tAKxMUwuIiJmnYk2zy0wxdUxT56SptUrTdUK6uGIAUmOqWd5E5cLgz3jazolNI7EIAhkAd5gabMnTbbq0tTqY2CNJG215jLvIacAhTMeNp23QpqYFEE1xcRNLTWnBdAQygEP40rUH1ypd2gJXTVxrZVKXqsFgIM+fP5dGo3Hnd77vS6PRkP39fcnlciIiUqlU1p1EIBWavlxcrARMBQ5/XftDA6nRxdT11pQ3WivdVbi4nIEJqq70YDCQp0+fTg1OgiCQw8ND6ff74yCmWq1Ku90mmIETTHy5aBpc66qHe/846ySo9funf7nU57RWkNq4FoCYouLbyvd9qdfr4nme7O3tTX3PycmJFIvFcRAjIlKv1+Xg4IBABvgaAUg8F1uZTKGCjGeq3NDimh4V33qe50mn0xERkV6vN/U9Z2dnUq/X73wuCAIZDAaSz+dTTyeQJte+XLTRFuRxve1gqtzQ4poea87G933xPO/O67lcTnq9HoEMrOfalwvS5+oy/NrSowXfEdNZkStBEMT+bm9vT96+fbvG1ADA6kw8obu6z5emncFNIThLj67SG+P8/Hzm72cFOgAgoq9ic5GmXau1jfuhHKVnY3J2NBpN/LyzsyM7OzsZpQbAummr2DTRtpAiQedmurq6kqurq/HPt+vtOFZd9WktL+fn5xMzmeI8fPhw4ufPPvtMPv/8c1NJA5AiKrbpXF1XxER6TLUOubjZqFYnJyfyxRdfLPw5K3IompI9rYspCAJ58ODB3GO8efNmYhtwWmMAe2ia8WGKixsjaspjbdslaLtWGh0fH8snn3wy/nk0Gt1phJjGikAml8uNp1pPUywW5x5jd3d3IpABkD5NT6Hanmy1pccEU+ekKSDC+iw75MOaO6lUKslwOJx4bTAYiIgw9RpQytRTqKZBpIhnavdrEwERwdDmUBfIBEEwteXl+PhYDg4OJAiC8ZiYVqs1XkgPgLtc3I3bRGuVqcraVKDo4nUiiNZPRSATBIGcnJyI7/vi+7602205Pz+XQqEgtVpNRN51L3W7XTk5OZH9/X0ZDodycHAgpVIp49QDiEMlEE/TOjK3W0RcYKo10MXgzDUqAplcLjd1t+vbPM9L9D4AOrhYCWga96MNgSuy4P6dBQAGMe4nnqZgTdsYGW3pcYmeUgfAOa7uBWSCptYqU4N0TQV5JtKjKagS0Zcel5CzAFLj6l5ArtEUVInoSo+mIJpuzencPjsAUEpTBYl4moIAFtWbTs8VAgALmApANO3w7OJ4HWwOAhkAWICpJ3QTQYipJ3RT0681tTKZGvejiab81YRABgAyoKmbQNvYCxMVtqZxNqZoCqI10ZkqAMBcpp7QtW2MqLXCXIWm4EFTEG2Ce6UFACxgIghxscI35d/+L7+fdRImuBY8aLIxd0GhUJDt7W05OjqSo6OjrJMDwFKmnqw1DfZ10Q//779e6nNprGmjjfaxNs1mU5rNplxfXyd6v/ul+WuvXr2S3d3drJMBpMpUxUYFGU/Tk7WmtLjK1FgbTcGD9vs0anAYjUZy7969ue/XfTYAFqJtrIMmBGfxTE2/NlVZm7hW2qaUb0I5ygo5C2AjmArOND1Zm6qstXWDmLhWps5JW0CEuwhkAIewt1H6ND1Za5ti7GJLnrYgD3fpuSMBrIy9jeIRnAFucu/bCgCm0BacmRgHQnAGEMgAmMLFgbHaptGa6IbRnN9Z03a9XbyntCCHANyhbayDiUpA23gSxHNxiwJt95RLCGQAqOdiJWCisnZ13SBNrRDa8gZ3bYVhGGadiDRFC+pcXl6yIB6QkLYv70effrnU524GMtq6GkwwkS8mj6Op3Ji63i7mjS2S1t+bm0MAYmn78jTRemEqIKFCiqfpHLUFoJryxjXkLAD1NFUCmrq5tM1acrGriwXx9NPz7QAAsJqLW2RoGzRsgqZA0QSdqQLgBNe+MEV0tYJoqvC1cbHsmeJauXH/in2tUCjI9vb2eFdNAOlz7QtTxFxFp6mi1ba1hYnuHG1lT9P1NiWtc2o2m9JsNuX6+jrR8fTmkGGvXr1i1hIANTRVtNq2tnCxO8fE9TYVOJgKONMqw1GDQzRraZ6NCWQArJ+mbhgX/XXtD7NOAtbIVOCguZVnGawjAwAZMPF0bWqNE21M5I22GVQuXu+0u8tYRwYAUmDqy9u1p2KTNOWNpm43bS2cWq6TjlQAgCU0jW0xtcaJplYHU0xdJ03npCVw0IZcAQBLmVq9Vtv6L5qCB02BK6YjkAGABWhr3neRieCB67Q5CGQAqKfpCZ11ZOxgKp9czBvXEMgASI2pytrF5n1N56St0jcxrVxToIh0ccUApEZTZY30mQoCfv/0L5f63M1yo23cD9JDIANAPRNP+jyhI0ua1sZxDQviAUiNpi/eTVtMLAumzsnEtdI2pdzEOWkrw2ljQTwAmXNxYKwpmtO2LE3dMC7mL6bjSgNQz0QFaWrxOFNcDM5M0TRTyFRwZqL8acoXTdy/I75WKBRke3t7vKsmgHi/+OWvlvrc7QXaNFXW2nZU1tR6YYq2YFETE+VvE4JYEZFmsynNZlOur68TvX8zckVEXr16xRgZOM9U4GCq0jdVWfMkagdN5cbUvUDZW7+owSEaIzPPxgQywCZw8SlfxM0N97SlxzWm7gUXW0E0tZSaoDNVADJlqotAU2Wt7UtYW3pMMFVuNHVRuVbpi7j3wKM3pwEszFTgYGozQs1f5stysWIzxVTXkonjmLoXXKv0TdJyL7h/ZwEbZBMqy6xRsdlB272gqXXStSBP15UGACSm5Yk4oqmy1jbY10SemzonbUHeqtw6GwDYIFqeiCOaKkgXB/tqu95aAlc9VwgANoi21hQTtLWCIF1ayqKOVCzA931ptVry4MEDefv2rTx48EBqtVrWyQKwITSNLzCVFlMBiLYWA9cQ4E1nVSATBIG0Wi1pNBrj13q9nlSrVWm1WhmmDMCm0PIUapK2AERTekytcm2Ci2XPBKty5eTkRKrV6sRrxWJR6vV6RikCsA4udsOYeLrWVOGLuNlioG1rC9yl9y6fwvd96fV6UqlUJl7f29vLKEUA1kFbhW2CpiBL08wcETPp0bSoHtKl505KoFAoSLValb29PSmVSiIiMhgMJJfLZZwyAFg/bQGIKSbSY6prx1RA5GKrohZbYRiGWSdiEQcHBzIYDKRUKkm1WpVutzsxZua2aNOpy8tLNo0ELEUlkC7yN56pvHn06ZdLHUdzq2JZGaKJAAAbg0lEQVTaktbf1pXCfr8vh4eHcnZ2JmdnZ9Lv9xN9bjQaTfy8s7MjOzs7aSQRgGEuVpiaggdtXXcmBtgyE8s+V1dXcnV1Nf75dr0dx7pvh3q9LuVyWcrlslSrVTk4OJBWq3Vn3MxtDx8+nPj5s88+k88//zzFlAJAPCrIeCYG2GrLXxcHQpt2cnIiX3zxxcKfsyqQqdfrsr+/Pw5aisXiOKB58uSJ5PP52M++efNmommK1hgAwDyujkPS6Pj4WD755JPxz6PR6E4jxDRW5ezZ2ZkMh8Pxz57nSb/fl4ODA3n+/PnMQGZ3d5cxMgDU4Ak9nqYZRwQg67PskA9rrpDv++J53tTfHR8fy6tXr9acIgBYnqYKUlPgIGJmxhGB4ubQcyfN4Xme+L4/9Xfn5+dSKBTWnCIA82ga0Ip42hZ9M1FuXCxD3E/TWXV21WpV6vX6xHRr3/el2+1Kp9PJMGUAptE24BJ2MFFuXKz0uZ+m03vFpqjVanJ2dibVanW8CN6DBw8IYgBgBS52w1Dpbw6rAhkRkVKpNF7VF4BumipIU0/oLj7pm0qbi3mjiab7SRNKD4DUaKqgTD2h86QfT1PeuFjpa7qfNCFXAACqmAhCtLUymTiOiRWPTaVFE52pAgDDXHxCd5WmClNTS56p2WWaWs5M0FNaACCGpum4pgIi156KtSF/NwdXDIB6mp4gTVV0ms7JRdry10QAbGrhQtdaJwlkAGwEF5/QNY3fEHGvgtTGxIrHIrrL9DK2wjAMs05Emkajkdy7d08uLy/ZawmwlImK9tGnXy51DM2r15o6J215Y4K2IM/FPE5b0vrbrbBshkKhINvb23J0dCRHR0dZJwfAAlx7ghRx85xcxHVav2azKc1mU66vrxO9nxYZABtBW9eSifRoa3XQRFsLiIt5nDZaZADgBm3ripgYjGrqnDa5soT9KL0AsABts2FcpKn1QlPgiukIZABsBE2Vo4iZGT7azskUE5W+qanKBCD66S7NAGCItgrJRDCh7Zw0MbUKLvQjkAEAqOLiejSmWohwF7OWAKjHDJ/pNKVFG23XW9ssKhswawmAM5jhM52mtGjj4vXGdFwhABvBxdYLF88JWBSlGYB6JsZMuDgw1sVzAhZFIAMAQMpcHMCsBYN9AaRG00BJF7thXDwnEV2DuzVx8ZxmYbAvgMxp6vqw9ct8FhfPScRMudFU9kxx8ZxMcPMuAOAUmuUBxCGQAZAaUwGIqy0PmM5EuSH43RwbM0bmBz/4gWxvb8vR0ZEcHR1lnSwAltq0cQrQY1PKXrPZlGazKdfX1/LTn/507hiZjQlkGOwLwARWaIXtbAmIGOwLAEjElootC7/45a+W+tx3vvVNwykxx7VBw+6XQgAwyNTYC03Bg2sVm0nsoq0fgQwALMBUIEHwAJhBIANgI2hqAdGGGT7xXv/Z06yTYJxr19v9OxQARF8LiKbKZBOCtWVpHuuyLNeut1tnAwCWcK0yAbLCnQRgI2hqATGF7jKAQAbAhnCx8tbWXQZk4b2sEwAAALAs9x5RAABwEF2J07l9dgAAKGAiCNHWlaglsCKQAYAMmKgEXBzA7CptQYgJWs6JQAYAMmCiEnC9ywCTCFyn25i7oFAoyPb2thwdHcnR0VHWyQFgKS3N6bCLiSBEWxlKK7BqNpvSbDbl+vo60fu3wjAMU0mJEkm3AQeAJB59+uVSn7vdnG4iICKogsuS1t+UZgCwlJYxCkCWCGQAYAGmmtNdDEJoIUIWKD0AsAAXK11TAYiLwRn0c++OBAALmGjZoXUIIJABgEyYaNnR1jrE9GBkQdddAABYOwIQ2IxABgBgBF1UyIKVgYzv+9JqteTBgwfy9u1bKRQKUiqVsk4WAFiJAAQ2sy6Q6fV60mg0pNPpSC6XE9/35fDwUIrFouRyuayTB8BxTDGORxcVspD4zvrbv/1byeVy8ujRoxSTM1+5XJaXL1+Ogxbf9+X8/DzTNAHYHKZaLzQFRKYCkE0I1qBP4lLX7XZla2tL/vRP/zTN9Mx0enoqnudJPp8fv1YsFuXi4iKzNAHAMjR15xCAwGYrld5PP/1UvvOd78QGNz/5yU/kww8/FM/zpNFoyB//8R+v8uek1WpNBDEAkJSpFhC6TwBdEgUyz549k263Kx9++OHE6/v7+7K1tRX7uY8//lgqlYr8xV/8hXz66acrBzK+70upVJJ2uy0iIkEQiIhIrVab+9nRaDTx887Ojuzs7KyUHgD2MNUCYqr1wsWASFN3GexzdXUlV1dX459v19txEpWejz/+WJ4+fSo/+tGP5KOPPpJvf/vbIiJyeXkpb9++nfqZn/zkJ3J5eSn1el1E3gU9q4iCFt/3pVqtiud5IvKuu6lcLkun05n5+YcPH078/Nlnn8nnn3++UpoAYFmaKm+2KIAGJycn8sUXXyz8ucR3kud58kd/9Efy6NEj+eijj8TzPPnxj38sh4eHsQkqFovjwcF7e3sLJ+6mmwN6oyBGRKRSqUi9XpfBYDCz2+nNmzcT24DTGgNsFhdbQEwhAIEGx8fH8sknn4x/Ho1GdxohplnokaDVakkul5Nnz55JEATSaDQkn8/L8fGxnJycjN/38uVLGQwGMhgMxq/9yZ/8ySJ/6o4oeCkUChOvR7OXer3ezEBmd3d3IpABsFk0tYC4imARq1h2yMfCd3aj0ZBGo3Hn9R/+8Ifywx/+UEREfvzjH0u1WpXf+73fWzhByxoOh2v7WwDWi7EX6WL6NWxmpNQ9ffpUnj59Ks+ePZPhcCjPnj1buQVmmnw+HzsmZ9UxOAD0ousjXQQgsJnR0vvxxx+bPNwd1Wr1zqBe3/dFRNiiAACWRIsXbLYVhmGYdSIWsb+/L51OZzweplwuy97enrRaranvH41Gcu/ePbm8vGSMDGApKtp0Pfr0y6U+R4sX0pS0/rbuLu/3+1Kv1yWXy0kQBFIoFBKtIwPAXgQkAOJY1yKzKFpkAGA2WrygkbMtMgAAswhIYLP3sk4AAADAsghkAACAtWhPBIAFMJ4E0IU7CwAWoGlxvl/88ldLfe473/qm4ZQA2SGQAQBLPfnzl0t9jvVf4BICGQBYgKl9ieiiikfeYBEbc9ULhYJsb2/L0dGRHB0dZZ0cAJYyVVma6KJ6/WdPjaRFG03dd1i/ZrMpzWZTrq+vE71/YwKZV69esSAeAKcw1gUuihocogXx5tmYQAYAYAdT3XfYDAQyALAAxm8AunBnAcACTI3foNUhHmNksAgCGQDIgIstNLRWIQvsfg0AC6Cyjvfo0y+X+tztlhTyGCLsfg0AqaCyTB95jEVQWgBsBJ7yATdxhwLYCAwgBdxEIAMAMIKZWMgCgQyAjeBiJautu4xuOGSBUgdgI2irZE0EIXSXAQQyAJAJghDADAIZALCUtu4ybV1d2AwbU3oKhYJsb2+Pd9UEgCyZCEK0BQC0MsGEZrMpzWZTrq+vE72flX0BAEawsi9MYmVfAEAiv/jlr5b63He+9c2Jn011ddGyg0UQyADAhnvy5y+X+tztwIEWEWSBUgcAUEXbIGboRiADAJYy1SWkDS07WASlBQAsZapLiBYQ2IxABgA2HC0gsBmlFwAs9frPnmadBCBzBDIAkAETa6VoH+sCrAOBDABkgLVSADPeyzoBAAAAy6JFBgAWYGr5fGYKAWYQyADAAkx1CTFTCDCDriUAAGAtHgkAYAF0CQG6bEwgUygUZHt7W46OjuTo6Cjr5ACwFF1CQLqazaY0m025vr5O9P6tMAzDlNOUqdFoJPfu3ZPLy0vZ3d3NOjkAACCBpPU3Y2QAAIC1CGQAAIC16OwFADjH1Ho/0I8rBgBwDltAbA66lgAAgLWsb5E5PDyUbrebdTIAAIqw3s/msDqQabfb0uv1sk4GAEB0jUthrMvmsPZKB0EgnU4n62QAAL7GuBRkwdoxMu12W6rVatbJAAAAGbKyRcb3ffE8L+tkAABuYFwKsmBlIHN2dia1Wk3Ozs6yTgoA4GumxqVoGmsD/ay76r1eT0qlUtbJAABnaAscGGuDRVgXyAwGAykWiwt/bjQaTfy8s7MjOzs7ppIFANYyFThoC4hgl6urK7m6uhr/fLvejmNV6Wm321KpVJb67MOHDyd+/uyzz+Tzzz83kCoAgIi5gIixNpvp5OREvvjii4U/Z00gEwSBiIjkcrmlPv/mzZuJbcBpjQGAd7QFDrTQbKbj42P55JNPxj+PRqM7jRDTWFNaer2e9Pv9iSnXvu+LiEi1WpVcLieNRiP287u7uxOBDADgHVOBg6aAiG4u+yw75GMrDMMwhfSsRbSWzKxTGI1Gcu/ePbm8vCSQAeAUKut4jz79cqnPMWBYj6T1t9WlOepuAoBNpG12D4EVsmBl6fF9XxqNxnifpcPDQymXy0sPBAYArE5TYKWpmwvpsjKQ8TxPWq1W1skAgExRWcejlWdzcKUBwFLaKmsCK2RB110AALCWtsAKm8Ha3a8BAAAInwFgwzHbCDajFALAhtM02whYFF1LAADAWrTIAMCGY7YRbEYgAwAbjrEusNnGdC0VCgV5/PixNJvNrJMCANgw//kf/stS/zZRs9mUx48fS6FQSPR+qzeNTIJNIwEAyzI1o4tNLBe3EZtGAgCQJmZ06UcgAwBAyhhQnR4CGQCAKpoW6DMVgDCgOj3kLABAFU3dOQQg+m3MrCUAAOAeQk0AgCqMJ8EiCGQAAKrQnYNF0LUEAACsRSADAACsRSADAACsRUckAMA5mtaiQbq4YgAA52haiwbpomsJAABYixYZAIBzWItmc2xMIFMoFGR7e1uOjo7k6Ogo6+QAAFLEWBd7NZtNaTabcn19nej9W2EYhimnKVOj0Uju3bsnl5eXsru7m3VyAABAAknrb8bIAAAAaxHIAAAAaxHIAAAAaxHIAAAAaxHIAAAAaxHIAAAAaxHIAAAAaxHIAAAAaxHIAAAAaxHIAAAAaxHIAAAAaxHIAAAAaxHIAAAAa21MIFMoFOTx48fSbDazTgoAAIjRbDbl8ePHUigUEr1/KwzDMOU0ZSrpNuAAAECPpPX3xrTIAAAA9xDIAAAAaxHIAAAAaxHIAAAAaxHIAAAAaxHIAAAAaxHIAAAAa30j6wQsqtfrSbfblSAIxPd9KZfLUqlUsk4WAADIgFWBzGAwkMFgII1GQ0REgiCQDz74QPr9vrRarYxTBwAA1s2qrqVWqyW1Wm38cy6Xk0ajIe12W3zfzzBlAAAgC1YFMi9evJDT09OJ1548eSIi77qcAADAZrEqkNnb25O3b99mnQwAAKCEVWNkhsPhndeiLqWoZQYAAGwOqwKZaVqtlhSLRcnn8zPfNxqNJn7e2dmRnZ2dNJMGAAASurq6kqurq/HPt+vtOFZ1Ld12dnYmvu9Lp9OZ+96HDx/KvXv3xv9OTk7WkEIAAJDEycnJRD398OHDRJ/bCsMwTDltqQiCQA4ODqTb7YrnebHvG41Gcu/ePXnz5o3s7u6OX6dFBgAAPaa1yDx8+FAuLy8n6u/brO1aKpfLc4OYm3Z3d2dmBAAAyM6yDQxWdi1Vq1VpNBoTQcxgMMgwRQAAIAvWBTLtdlvK5fLE4F7f91kQDwCADWRV11Kv15NOpyOHh4cTLTDdbne8bQEAANgcVg32vX//vgRBMPV3cacRDfadN1gIAADokbT+tqpF5uLiIuskAAAARawbIwMAABAhkAEAANYikAEAANYikAEAANYikAEAANYikAEAANYikAEAANbamECmUCjI48ePpdlsZp0UAAAQo9lsyuPHj6VQKCR6v1Ur+y6DlX0BALBP0vp7Y1pkAACAewhkAACAtQhkAACAtQhkAACAtQhkAACAtQhkAACAtQhkAACAtQhkAACAtQhkAACAtQhkAACAtQhkAACAtQhkAACAtQhkAACAtTYmkCkUCvL48WNpNptZJwUAAMRoNpvy+PFjKRQKid6/FYZhmHKaMpV0G3AAAKBH0vp7Y1pkAACAewhkAACAtQhkAACAtQhkAACAtQhkAACAtQhkAACAtQhkAACAtQhkAACAtQhkAACAtQhkAACAtQhkAACAtQhkAACAtQhkAACAtQhkAACAtTYmkCkUCvL48WNpNptZJwUAAMRoNpvy+PFjKRQKid6/FYZhmHKaMjUajeTevXtyeXkpu7u7WScHAAAkkLT+3pgWGQAA4B4CGQAAYC0CGQAAYC0CGQAAYC0CGQAAYC0CGQAAYK1vZJ2ARfm+L41GQ/b39yWXy4mISKVSyThVAAAgC1YFMkEQyOHhofT7/XEQU61Wpd1uE8wAALCBrOpaOjk5kWKxOA5iRETq9brU6/XYz1xdXU38D7Ourq7k888/J39TRB6ni/xNH3mcrk3PX6tW9t3f35d6vX6n9WVra0v6/b7k8/k7n/n5z38uDx8+lDdv3sj3v//9dSV1Y7BycvrI43SRv+kjj9Plav46ubKv7/vied6d13O5nPR6vVT/tqk9mjQdR9u+Uy7mjaY81nZO2o5jgqZz0pQWk1zMG015bOU5hZa4uLgIRSTsdrt3fud5Xlir1aZ+7s2bN6GIhG/evFnp7//u7/7uSp/XeBwTx7i8vAxFJLy8vFSRHlPH0ZQWU3ms6Zw0HYcynP5xtOWxprwxcRxt+WvqOEnPy5rBvufn5zN/HwTB1NfDr3vO/v7v/37i9Z2dHdnZ2Un896+vr2U0GiV+vw3HMXGM6PNazsnUcTSlxVQeazonTcehDKd/HG15rClvTBxHW/4ue5yrq6uJcT5fffWViPy2Ho9jzRgZ3/dlf39fut2uFIvFid/t7+9LsViUVqsV+zkAAGCfeWNcrWmRiUxreTk/P5+YyXTTo0ePZDgcyvvvvy9bW1vj1xdtkQEAAOm53SIThqH8+te/lu9973szP2dNILO3tyci07uYgiCQBw8eTP3ce++9N3WAMAAAsJ81s5ZyuZx4nhc7FuZ2dxMAAHCfNYGMiEipVJLhcDjx2mAwEBGZWEPm8PDwzmd935dqtSqnp6fSbrel3W6nm1jHTcvjg4MD6fV6EgSBBEEg7XZbTk9PM0gdAGBTWDPYV+RdF9LBwcGdLQoODw+lVCqJiEi73ZZqtToxyjnucwcHB2xtsIRpeSwiE2OQRN61knW73XUmzVoHBwfSaDTkyZMnIiLy4sULCYJAarXa+D3sM7aaJHmc5D2I5/u+tFotefDggbx9+1YKhcL4uzn6PWV4NfPyeCPL8MoTvddsOByGtVotbLVa4/8jFxcXYbFYDG+fVq1WCyuVyp3j5HK5taTZJXF5HIZhWKlUwlarFTYajbDf72eQOnuJyMS/YrE48fuLi4vQ87zw4uJi/FqU30hmXh4nfQ+m63a7YbFYHJfR4XA4UWYpw6ubl8dhuJll2KoWmXlOT0/F8zwpl8sTrQXLbG2A6eLyOPqd01F/iqIWwiAIpFgs3imT9XpdgiCYWGLA9305ODiQi4uLdSfXSvPyOOl7MN39+/fl5cuX4zzr9XpSLpfl7/7u7ySXy1GGDZiXxyKbWYatmbU0T9z2BbN+F21tsAkX2oRZeYzV7O/vz2xiPzs7u7M5ajT4fTAYUIYTmJfHSd+Du6IHnJvlsFgsTgQolOHVJMljkc0sw1YN9p3l7Oxsop8wEjfLSeTdlO63b9+mmSynxOVx5O3bt+PB1Kenpwz0NSjLfcaAeVqt1tyHHMrwapLk8aZyokWm1+vFVrDLbm2ASbPyOBIN5ItUq1Wp1+sTr2G6KAjM5XLjMhl10xGMmzErjxd5D+7yfV9KpdJ4Nihl2Lx5eRzZxDLsRIvMYDAgUk1ZkjzudDoTP5fLZTk9PSVYTMD3fanValKpVKRWq8lwOBw3wxOMmzErjxd5DyZF5c/3fSkWi+O8E3n3HSBCGV5VkjyObGIZtj6QabfbifoDF93aAL+VNI9vi6b/0Ww8H0Fg+pLkMddhcTeDlJsPO5VKRc7OzsZrfWF5i+TxJpZhqwOZ6MLMCkaW3doA7yTJY5F33UhxAYvv+8bT5bppQSDBuFlJAm2C8fmiirVQKEy8HpVLyvDqFsnj2zahDFs9RqbX60m/35dqtTp+Lao0q9Wq5HI5aTQabG2wgqR53G635eDgYOKzUfDIbITZqtWqlMvlqWUxakoWIRhfxbw8TvoeLG44HPJAmbJoxftNLcNWBzKlUunOANR2uy29Xm9irYKkWxvgrqR5HPXJ3tTr9SSXy42fCDDdvCCQfcZWlyTQJhhfXj6fjx2wG63iSxlezbw8FtncMmx119I0026U4+Pj8R5AkVardacvEclMy+PDw0M5OzubeE+j0ZBnz57RbDxHkiCQYHw1SfKYYHx51Wr1zliYqAUgehCiDK8mSR5vbBnOdmFhc4bDYVipVELP88bLMt9c+nrW1gZIZl4ed7vd8XYQpVIp7Ha7GabWHt1uN+x0OuOfo6Xcp712e3n3m+9BvCR5nOQ9iOd53sTWJKVSaWJrGMrw6ubl8aaWYae2KABs1ev1pNvtShAEcn5+LtVq9U5ze7RZ3P7+vgyHw41cwXMVSfI4yXswXRAEUq/Xx+uX7O/v31m/hDK8miR5vIllmEAGAABYy7kxMgAAYHMQyAAAAGsRyAAAAGsRyAAAAGsRyAAAAGsRyABw2uHhYdZJAJAiAhkAzjo9PXV6szwABDIAlKvX63JwcCBbW1uyv78v5XJ5YjuMOIPBYLx0u8sb5gGbjkAGgGqNRmO8+3q325VOp3NnI9M4nudJLpcjkAEcRiADQL1utyue54nneYnePxgMxhsRep5HIAM47BtZJwAA5un1eon3i/F9X3q9nrx+/Xr82u1dlwG4g0AGgGpBEEgQBIlnHwVBMLGR3nA4pEUGcBhdSwBUi2YdJWmRabfb4y6lyP7+PoEM4DACGQCqdbtdyeVyM8fH+L4v5XJZGo3GRNDS6/Wk1WrJYDCQ09PTdSQXwJpthWEYZp0IAIizv78v+XxeOp1O1kkBoBAtMgDUCoJAfN+XQqGQdVIAKEUgA0CtaOZR3PiYdru9zuQAUIhABoBa3W5XROTOAN4I06oBEMgAUGvW+jHtdpsuJwAEMgB0CoJgYoXem9rtttTr9cRbFQBwFwviAVCnWq2O14/p9XrjvZZ835fXr19LEARSqVTufK5er0u73ZaXL1+K53lyfn6eeFsDAHZi+jUAJ/i+L77vy5MnT+TFixfy5MmT2LE1ANxBIAMAAKzFGBkATgqCIOskAFgDAhkAThgMBtJut8X3fQmCQE5OTrJOEoA1oGsJgPWi/ZX29vbk6dOnIiLS6XQY6AtsAAIZAABgLbqWAACAtQhkAACAtQhkAACAtQhkAACAtQhkAACAtQhkAACAtQhkAACAtQhkAACAtQhkAACAtQhkAACAtQhkAACAtQhkAACAtf5/pjW0KXcFhg8AAAAASUVORK5CYII=",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x151710fd0>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = subplots()\n",
    "ax[:plot](ent_spec[:, 1], ent_spec[:, 2], ls=\"none\", marker=\"_\", ms=\"10\", mew=\"1.5\")\n",
    "ax[:set_xlabel](\"\\$L_z^A\\$\")\n",
    "ax[:set_ylabel](\"\\$\\\\xi\\$\")\n",
    "ax[:set_title](\"\\$N=$(N), N_\\\\phi=$(Nphi); N_\\\\mathrm{orb}^A=$(length(subsystemA)), N_e^A=$(NA) : P[0|0]\\$\")\n",
    "ax[:set_xlim](40, 68)\n",
    "ax[:set_ylim](0, 12);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculating entanglement spectrum for NA = 8, LzA = [50.0, 51.0, 52.0, 53.0, 54.0, 55.0, 56.0, 57.0, 58.0, 59.0, 60.0, 61.0, 62.0, 63.0, 64.0, 65.0, 66.0, 67.0] ...\n",
      "elapsed time: 446.294034166 seconds\n",
      "\n"
     ]
    }
   ],
   "source": [
    "subsystemA = collect(15:get_Norb(setup))\n",
    "NA = 8\n",
    "LzAvec = collect(50.:67.)\n",
    "\n",
    "ent_spec = entanglement_spectrum(psi, setup.states, subsystemA, NA, LzAvec);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjIAAAHUCAYAAAAgOcJbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3d9rY2l64PHH7e7xMExslWsSmGEM1UfkImZZgmzB3mxYKHn2IrCwidz9F7S04L1qgtSGIdOzhPHKzK1IpPoLquTOXV8MUsGSvVhIlUQulmVvdDJQww4kU/axewjjTJyzFzVHY9mSfCS9R+d53/P9QNFt+ej4Pe/58T7n/bkWhmEoAAAAFnov7QQAAAAsikAGAABYi0AGAABYi0AGAABYi0AGQOqq1ar0er20k4EHcJ6gEYEMgFT1ej1pt9sSBEHaScEMnKfk9Ho9efTokTx69EgODg5W/vdPT09Hf79ara787y+LQAZAqrrdroiI+L6fckowC+cpWR999JFcXFyM8nmVarWaXFxcSKfTWfnfNuH9tBMA+wwGAzk5OXnwoo+2i3ieJ41GI+nkOaPdbku32x0VHB9//LHUarWJ2w4GA2m1WiIicn5+LsViceq2mpyenkq1WpWzszMZDodpJ2cpLt8X2s9Tr9cby/dcLifHx8eSy+VSTJUZ0bHlcjkJgkDy+fzEe9v3fWk0GqPnwCRxr1HrhFiJcrkcFgqFUERCEQkvLi6mbpfL5UIRCT3PCyuVyopTOt3FxUXYarVCEQlLpdLMbVutVlgoFMLhcDj6rNPphLVazUhaXMjPWWq1WtjpdEY/D4fDMJfLhZ7n3du21WqFnueN5UG5XA7L5bKx9CSR38PhMGw0GmEYhmGpVDKa3lVy/b7Qfp5arVbYarXGPut2u2GhUFhpOsrlcuh53ijvo7wql8thqVQKC4VCWCqVwn6/f++73W534jnodDr3jqPRaEw8B/1+f+r1F/canZYO7QhkVii60ERk5oOr0+moe1hEN2Or1QpzudzMm6Hf74e5XO7eQ9TzPKMPF5vzc5Z+vz/xeLrdbigiYw+ai4uLUETGgp7bn3e7XaPpMpnft48jKoBtk4X7QvN5uri4mJrnjUZjFICtSr/fv3eP3hYFE3fv12kBRC6Xm3gP53K5e/uYFsjMc40SyOBBjUZj9Gady+WmbtdqtSZG7Vo8dDNMe8MzXSPiSn7eValUpr5JR2/Rt7ed9ubted6DNQTzMJnfrVZr7AFdq9Vm7tMGLt4X2s/TrIJ3Vg1FUqYFKrfdvYfDcPJxdDqde9tFJtW4xjleVwMZOvuu0HA4FM/z5KOPPpIgCKYOYxwOh1IoFFacOjN6vZ74vi+Hh4f3fjccDme2387L1fx8/fq1fPjhhxNHh3ieJyK/7XDZ6/Ukl8tN7AtQKBSMDpU1ld9BEEir1ZJutyv1el3q9boMBgOnR8PYeF/YcJ5mHd/5+fnK+8hEHXVLpdLUbba3t0Xk4U7T3W53avo9z2MY/C0EMimo1+siIlM7+Gl6UMwr6kQ260Y2zbX83N7eliAI5Pz8/MFtZz0M4z4w57VsftfrdXn58qU0Go3Rv2jIp23nKi4b7wsbzlOhUBDf92Vvb+9emlqt1sqHEvd6PfE8b2YAFd2P0f350L4myefzEgSBDAaDxRPrEAKZFRkMBrK3tyci76Lp6G357s0X3ZS2uv2WcHp6Kqenp6O3OZNczs9OpyP9fn/iQyx6CE57wN0WPUxvP+yq1ark8/m502Qqv9vtthwcHNx70N+taXKNbfeFLefJ8zwpl8syGAzkww8/HOVzu90Wz/PmChwXvTcivu9LEAQz/+bZ2ZmIiJTL5Qdri+K8yMTZJgsIZFak1+uNXeDHx8ciImPDMCdtZxvf9yWXy0m73ZZKpSK1Wk0ajYbk8/nRW4QJLudnLpebWOUfFUi3h17OahqYVNhETRzzFkTL5nfUrFKtVu/Nk3F2djYq0Ov1upNV5rbcFzaep06nI+VyWYIgkIODA9nb25Pt7e25h7Qvem/c/r6IzAwQ6/W6FAqFWMOfgyCYWmsTfa6lZix1aXfSyYpJHahE5F7HuUWGYUYd8Bb5t8iollkdxqJjujscMgzfdWo0Ncw0yfx8yKrzO1IoFO4Ns67VajM7+4rI2LkYDocLdXxOM78XxX0x/rdvW9V5WtU5uLi4GHWild8MKb89xD2ORe+NSLlcDkVk4t/t9/ujTt2T7tVJnWxlxlDpqCPw7espy519CWRWZNLFEY04uX3D2nARPfTAnlawRjf6tBE58zCZn7NGGGjRaDQmPiSnDb8eDoejIbgmjs+l6zcprt0Xtuh2u6HneWG32x3dJ1HwNm8ws4wokKrVavf+RSPJZh0DgcziaFpagWnt0Xc74/m+v1QbrSaT2n9Nta+bzk/t7cy+78vJyYl0u917fWNyuZx0u105OTkZVTMHQSBnZ2ejbZcduZHF6zcpNt0XNoiawrrdrpRKJanVaqP+ZUEQTBwlllQ6ov4xtztHR/9qtVqsfm1xRc+shzoMZwVLFKzAtPbou53xBoOBdf057ppVaD5+/FhE3g0vXmY4dJbyU0Tk8PBQOp3O1GMplUrS6XSk3W6PPqvVaqOH+P7+/lJ/P438joIyF6aYF+G+SEq1WpXj4+OxIKFQKMhwOJSDgwPp9XoyGAwSn34h6h9jMnDK5XIPvmS5cn8sixqZFZg2AkVERsMD2+22vHr1yqr5TibZ39+f2gHt7du3IhJvxM0sJvPT932jb0qmHR4eyvHx8YMFj+d5UqvVRv9Efntsyz7sVn399no9efr0qTx9+nTpfWlh231hi16vJ+VyeeLvonlYVjHCKs78MfOKapUmida70vzsWiUCmRWYVZBUKhURkbGmgXnV6/XREuzz/jM98iBagn7SsUSfLXvzmcjPdrst7XZbGo2G+L4v9Xo99gNvVfldr9fl4ODg3oO63W7HSutgMDAyvDfp6/euqHp+WdwXv5XEeYpjFedgVvPK/v7+SmotookpTQYWH3/88dT7PAgC43/Paml30nHdcDicOFLhtqgz3qrXBVnUrA5jw+Fw6ho/hUJh4jot83RyNJGfpVJplL5oXxcXF6rWkOl0OlPTf3tq8m63G+ZyuXujLaJ1Ve6ad2RGWtdvGov+Lcv2+8JGt+/lSe6O8pvl4uJioVFL0bldZj23SZ1so/1OStOkZS3o7IvENBqNB6PmqNrXlnbrWbPOep4nlUrl3hv1YDCQwWAgz549G/t8b29PHj16FLs2ZNn8rFarcnBwcO93uVxOzUyZvu/LJ598Is+fP5e9vb2xf/l8fiyNk/ItCAKp1+v38lpERvNsxH0rd/H6TYrN94XIu2vp9PRUzs7ORv+0uz3b8F3tdlvq9XrsGpnoHpu3KSrKp6jWzZRoQr/nz5+PfT4YDMT3/YVqLePOGG6dtCMpV1UqldE8Hrlc7sFofdWLm82rVquFpVLp3lwNpVJp4vDeUqkUViqVsN/vh61WazQ88q5yuRx6nhfrbdJEft6+5Pv9/tjbjhgaArusQqEwyuNJ/+4eW7lcDjudzlheTxtyXalUYs3RYSq/Ly4uwlqtFnY6nbDT6YS1Wm2Ux1GtS3QNVSqV0dtg9Lto0cJWq6VqjpqIK/dFNFz/9s8a83uSTqczulb6/f6oVmHemqcoL+Oq1Wpj96rneWG5XF5oyPe0mpCopji6Dvr9flgoFCZeM9NqZOa5Rm2tkSGQQWKih3Wn03kwQHjogW1CtCLwpL/ZarXUB5OzRE1RrVbrwbyOHvarUCgUxoLF4XA4Vlh0Op3R5GX9fn9U+HS73XuBZafTsfocRbTdF2H420n5Wq3WKG22ieaRiZOv2jwUQMQ5NhOrfdsayDD8GokpFAqxRkcMBoOVdFrzPG+0IOPt6ubBYCCdTifWtOFaTRu5Mcn5+flK8vvs7EyCIBi7BqJRVNFU/dEQ0yg9k7aNlMtlOTw8XMlw2iRpuy+CIFi4qUKTUqnkbPOmy8dmAn1kkLrnz5+v7CaNJo87OzuTfr8v7XZbXr9+PRqqmQWrKiBfvXo1MU+3t7dHw0dF5h+to2XBwqSt6r6YNTIGsAGBDFJ1dnY2tbNeEjzPG3WMrNfrUqlURkNTs2CVhVOxWJxYSC47s2wWhpyu8r4oFAqjzu4R3/fVLAoJPIRABqkql8upFEyvX7/ORIF414sXL8ZWz05SuVyW7e3tiaOs4gSP0bTvkXa7LaVSyepmpbhWfV90Op1RTeXZ2Zn4vj9XcyWQJvrIABmy6tqnfr8vJycn8vr1axF5NyNpv98XkXeTiDUajdHEfdVqdazwbrVao1qB8/NzGQ6HoxlUYRZ9MNL34sULefHihezv7y90nS8zQd7p6amcnJyIiMhHH3200D7StBaGYZh2IgAAABahrmlp1rTqvV5v9OZ2cHAwtkgeAADIHlVNS4PBQJ4+fTqx+juaATMaIhgEgXz44YfS7/el1WqtOqkAAEABFTUyvu/L4eGhPH/+fOoCYK1Wa6yTYi6Xk0ajEXsBPQAA4B4VgYznedLpdKTRaEydy+PFixdyeno69tn+/r6ICMMEAQDIKBWBTBzb29vy9u3btJMBAAAUUdVHZpbbM4FGoialqGYGAABkizU1MpO0Wq3MTJAFAADus6ZG5q5o9slocq1p/vVf/1V++tOfygcffCBra2ujzzc2NmRjYyPpZAIAgBiur6/l+vp69HMYhvLrX/9anjx5Iu+9N73excpAJggCqdfrsRb6++lPf7rUui4AACA9b968ke9+97tTf29lIHN4eCjdbjfWdMwffPCBiIj87d/+rXz7298efT5vjUyxWJRXr17Nn1jF+zGxj6urK9nZ2ZE3b97I5uZm6ukxtR9NaTGVx2kc07/5wU8W+hv/+4f/MZH0TDJv/po6JhP7+ad//pep2/3RH/2R/M3f/M3E333ja+OP/qTPUxrPCVPHZCqPTexH07Vncj+RuzUyX331lezu7srv/M7vzNyfdYFMtVodrV4cGQwGU/vJRM1J3/72t2dGdA9ZX19f+gbUth9TaRER2dzcVJMeTXmjKY/TOKb3Nr6x0N+YJ52mjitu/v7fxn+e+ru9vb2pzd13CzYTefP+jMJxVr4kkZa428f5zqxC/72vfV3e//rk9N4+LlPH9G8/+3L6xv+pIf/ux/9r4q9++t//eOxnE+fK1DFp2880V1dXIiJj3UImsSqQabfbcnh4OBa0+L4vvu8n3uH36OjIuf2YSospLuaNpjxO45j+z3+b/sbearWlWl1+EctV5/HdIOC2//pfKjN/f5uJvNn98xlvxH/8o6m/v1vIamPiuFZx7c0jyWO6uvpKvvOdb8v/+38/l83N2bUXs/YjMl/eaMljdYtG5vN5KZVK95YdiFbKPTg4GPu82+1Ko9GYGsj87Gc/G1VpLlMjg8murq5ka2tLLi8vjdU8YBx5PPsNfZbbQcW0fTxUCMQNTNLwZFZtwQx3AxkT+TtrP/PmsanjMsFU3pg4Jm3XsKm8mSbus09FIBMEgZycnIjv+3J2diYiIuVyWYrF4mhZgkePHkkQBBO/P+sQ/vEf/1F+7/d+T/7hH/5Bfvd3f9d84jPu+vpaTk5O5Pj4mFFgCSGPzRQC2gp9E5IOQObdj6k8djGQMbEfbddw0ucpbiCj4lUjWjdplouLi4X2HT34s1oAJG1jY0M+//zztJPhtHnyWFMh66qZTQQzJFGYmDpvpo7JRabyRtM95tr51pOzAJbm2gMqMqstfpX7MMnVc2WCiXPlYlCv7RrWQu8ZA4Df0FS4aCpMTBXWmo5JxMz5NhUoasobTfeBiJ680ZUrAJai5cFimokC28UaEE3NXJjNxTzWkjYdqQCgqs+ENpqCEE1pMWX/L14u9L2kjsnEvWAqqDd1vjVdN6+//9T4PtPk5lMPsJCmBx2S52K/H1O4F5KlLXBdFoEMAPU0Fdim0mKi9sxUDZxrb+gi+vrIuNiBWUt6CGQAJTQV1tqYePBpCkC0+dY3v552EtQydb5d7MCspebMvTsSsJSLBaQmLuavljfiiIs1O6ZoOleu3QtuHc0MxWJR1tfX5ejoSNX6N4BGmh66rnJxJJapvhcu9h8yca60HVNSms2mNJtNubm5ibV9Zp46r169yuw6NcC8tBWQLgZWmvJYW/5q6j+kiYvHNElU4RAtUfCQbOQKAKtpKvQ10dbXQVONgbbgTBPXJlJ0/4wBmJuWB5RJmhb/EzGTx1kodBflYvBr6tpzbf0oHakAoIqWB1TERKGvbWIzE3msrdaBmYan03QNu0bvWQeA39BcQKXJxSYhEV0FtqmgStM1rO18L0tPzgKABUwVAppqHUztU1MBqa3/kAnMhTSZW0cDQBVNhbUpptKmqYDUVutgYj4azdfQokwdk2v3pc5UAXCCi30mNKXFxVoHEV1rAblYA6ftfC+LQAaAepomE3O1X4oJmgprU1ysgTNFy/nWe/UAsJ6mwlpbYalp7R1TSwtoCvK0FLIaaQvql+X+GQOQGhc7kbpIU1OOKdqaNTUFZ64Fa24dDQBVePBmi7aaHRM0zRukLTjTQmeqADhBU4Gk6c1axEx6TKXF1H401exQizedpvvSBAIZwCHapuHXRNObtYiuwkTbedO0fIOpWiYTtAVnWtKj6+pNULFYlPX19dGqmoCLtE3Dr+VB5ypNwZCIrvNtKhjXVMukrc9ZUgFws9mUZrMpNzc3sbbPTCDz6tUr2dzcTDsZQKZoetPXVMiK6EuPJiYCNG1Bniaa7stJogqHq6sr2draenB73UcDYC7a+kyYYuLtWtusqCbSo6nZQ8TN4EFTnyhMRg4BDjH10NP28NRUQGoaOWKq2UNbIaspQHOxT5RrdD2tACAjNBVs2vpEaeqXoom2gFMLt48OgBO0NXVpoanmQhtt14yJ9GgKfjUhkAGgnotvlCYKNm01F5oKWm3XjLb0uIScBZAJ2ibEo2DDvLTVMmnBnQQgE7RNiKdpZl9TtKXHBE39Ugh+JyNXAGAOpgo2E4GVi7PXiugKiDQ1l2EyAhkAmWCqcHSxYDPV10ZTrZemmhQkizMGIDGaChNtBZS2WhDXaBtSjuTourMBOEXT5HGmuDhXiqszQpugLQDGfZwhAOppas6hYJtOU96wtEB2kNMAEuPiG7opmvLGxZozUzQF0S7mrwluHx2AVJl6gGoq9E1xsXDRFBBpCkBMcfGYTHDvTpqiWCzK+vr6aHlwAPZwsdDXRFugqKnA1pY3WdBsNqXZbMrNzU2s7dfCMAwTTlOqrq6uZGtrSy4vL2VzczPt5ACwnKnqfRP70dbU8Itf/mqh733rm18f+/nJZ18utJ/bgYy2vHHxfCctbvlt59EBQEpM1RaY2I+mmgsRcyOxTNSCaCu8NU2A6FpApDNVAAAgEdoC4GURyADAHEz1mdA0IZ62BTVNFLTaah3oa5McAhkA6mkqlEztU9OEeK69oYvoOyZNzTKuBVV6chaAczQtsOgiTbU6IiwLYAtNQZUJ6o5mMBjI8+fPpdFo3Pud7/vSaDQkn89LLpcTEZFKpbLqJAKIiQBkOhOFtabOtSaZKGi1HROSoyqQGQwG8vTp04nBSRAEcnBwIP1+fxTEVKtVabfbBDOA41wslFx7KxbRdZ60jfDR1DzqGhU55Pu+1Ot18TxPtre3J25zcnIipVJpFMSIiNTrddnb2yOQAZQyVbDxMJ9MU+dabbQ1a7qYx1qoeDp4niedTkdERHq93sRtzs7OpF6v3/teEAQyGAykUCgknk4A89EUgPBGnDyWKEAarLlDfd8Xz/PufZ7L5aTX6xHIAJjJxYLNxWMS0XVcpmq9NDW7ucaKQCYIgqm/297elrdv364wNQCwPE01RBSyyaNmLzlW5Oz5+fnM388KdADMT1Mha4q2wtpErYO2Pkia8pj+Q9mh96lj2NXV1djPGxsbsrGxkVJqAN14eNtBW+CoLT2wy/X1tVxfX49+vltuT2PVVTep5uX8/HxsJNM0Ozs7Yz//4Ac/kM8//9xU0gAop6kjqoiu2gttTEz0RzBun5OTE/nhD3849/esCGSiIdmTmpiCIJDHjx8/uI83b96MLQNObQwwHYXsdKYKSGovptO0fANW5/j4WD799NPRz1dXV/cqISax4k7K5XKjodaTlEqlB/exubk5FsgAmM7FQpbgLFtcHG3kYt+12xbt8mHH0YlIuVyW4XA49tlgMBARYeg1gAdp69CqqVDSNnutpuBBUxBAc9lkes7QbwRBMLHm5fj4WPb29iQIglGfmFarNZpIDwBWwVTBpqlQ0jZ7rYk81tYnSlPg6hoVORQEgZycnIjv++L7vrTbbTk/P5disSi1Wk1E3jUvdbtdOTk5kXw+L8PhUPb29qRcLqecegBIB4Vj8jQFeZpqqjRZC8MwTDsRSbq6upKtrS25vLykjwwANUwEIU8++3KhfSRV6/CLX/5qof1865tfN54eU8dkKo9N7SdL4pbfhOUA1HOx5kFT2kylxdRoI03Nbi52GnaNnjsJAKbQVLCZYiI4c7HjsSma+uuY3A/uI2cBIAUmCloXOx4D8yKQAaCeiZoHF2sdXGXifJuYHViEUUs2IIcAqKdpOK4pmvpMaEqLKdr662i7/lxCIAMgMbyFTmfiGE3lr7bJAin0MQ+GXwNIjKYhp9qCKheHX5ti4ri0NQlpymNNaZmF4dcAcIu2Wh5NtQ7aZsE1Qdv51pQeTdeeCXpyFoBzXOx7gelMFZCarhvXCn0XZSaQKRaLsr6+LkdHR3J0dJR2coBM0PQWqo2mwlpTWkR09R8yxdSsxyZoO993NZtNaTabcnNzE2t7+sgAgKU0FY4iuoIHbf2HNPUXswV9ZADAcaaGGJtiKiDRFBBRq6gfZwgAYISpAMTFlaJNTdCH+whkAGAOmmoLtBXWmjrGaqtJSao5DwQyADAXCuvkmai9cHH+F0xGTgMAjDDVfGKi7w9LCyRPS5BHIAMAc9DWnKOJts7HSJaWII9ABgDmQJNB8ggWMQ/uSACAEZoCEE1pEdHTDGOSljzWm0MAgJXQtoq2iSYLVvSeTtv5XpaOVADADC6+zWri4qKRpmhO26I0BVUmuHeGADhH04PXxcLaFE2LRmobfq2lGcZF7t9ZAGCQpqDKFG2FrImgT9vwa00LYWo738sikAGgnmsPXldxnpKlKajSxK2jmaFYLMr6+rocHR3J0dFR2skBMAdND14XC2tttUyamu9cPN/aNZtNaTabcnNzE2v7tTAMw4TTlKq4y4ADgG1MFfhPPvtyof3cDWQ07UdTMGTKL375q4W+Z+s6T3HLb71nDAAwk7aaFE1MBSSaAiJmTp6MQAYAYISLzTAEi/oRyACApUwFDqb2o7lZxgUuBoomcNUBwBw0NTWY2icBSPJMXDecp8nIFQCYA00NWATXTXIIZAAAzmHyuOxg+DUAzEFT05KrTOSxqaHgpnDdzI/h1wBwi7YVf10s2LTlsQkuHpNryFkAmaCtj4K29Jig6Zj+Z+0/GNmPpmPCZAQyAGApbSs8a/LvT//HQt/LQgDi2vnWmSoAMExbp00T6dG2wrMpmgpabdeNCdrO97IIZABkgra3SW3pMcFUoW+ioH39/adG0kKfKP3IIQBIgYmCjZl9pzO1UKKpAERTLYhrtUx6rrqEFYtFWV9fl6OjIzk6Oko7OQBWTNsbsYmCzdWZfU3UprgYgGRlBFWz2ZRmsyk3Nzexttd9NAa9evWKeWSADNNUIGmjLcgzscqztvOtqU+UdlGFQzSPzEMyE8gAgCYmCjYXax20cbHZzTXM7AsgE7TVOphgavZaF2fB1Xa+XTympDGzLwDcwuiT6bR1/tScV4vS1CfKFC33gq5cAQDlNDXDaAtANDF1nrQU1hppuRfcz2kAcJSpwlJLgaSRqbwh6EwOgQwAzMFUgcSbfrJMTYhniovnTUtwRmdfAEiBpg62Lq7ZZCp/NR1T1jjb2df3fWm1WvL48WN5+/atPH78WGq1WtrJAqAcBdJ02pqoOFeYh1VnPQgCabVa0mg0Rp/1ej2pVqvSarVSTBkA7bT1A9FSLa+RiXOlad0nEYKzJFmVQycnJ1KtVsc+K5VKUq/XU0oRACzGxQJKU3CmLX+1BdIu0XWmH+D7vvR6PalUKmOfb29vp5QiALbQVMi6ylTwoOlcaUoLJrMqkCkWi1KtVmV7e1vK5bKIiAwGA8nlcimnDIB22t7QMZ2L54qAKDnWjVra29uTwWAg5XJZqtWqdLvdsT4zdzFqCQCwKE2jy7LG2VFL/X5fDg4O5OzsTM7OzqTf78f63tXV1djPGxsbsrGxkUQSASCTNA3j1ta5Vlt6NLq+vpbr6+vRz3fL7Wmsq5Gp1+uSz+dFREYdf1ut1r1+M5Fpy4D/4Ac/kM8//zyxdAJA1mhaxFLbPDIm0uN6MPT555/LD3/4w3ufO1UjEwUxUdBSKpXk8PBQqtWq7O/vS6FQmPrdN2/ejGUEtTGAPVx/gEMvTdeQ6yOfjo+P5dNPPx39fHV1JTs7Ow9+z6oamXw+L8Ph8N7ne3t7UiqVJvaVoY8MYD/6KdiBpqXpfvHLXy30vW998+uj/8/afeBcHxnf98XzvIm/Oz4+llevXq04RQCwOG0FrQma07YoU+dp/y9eLrSfJCb5c401V53neeL7/sTfnZ+fS7FYXHGKAKyKiw9wTc0E2oIqE3ljKn81nSdMZk0gI/Kuc2+9Xh9rQvJ9X7rdrnQ6nRRTBiBJLr7pa0JhnTwTwTjnaTKrng61Wk3Ozs6kWq2OJsF7/PgxQQyAlTFVe+FiLZMpJvLGVP6a2g/BeHKsy9lyuTya1RcAVs3UW7Gmgk1bUGUib0zlr6bzhMk4QwCQcS4W1tr6/WjiWt5YNfx6EQy/BmCSa4WAq1wcqqxpcr5VcG74NQBoQECSPILFyVw/vkVRIwMAMEJTjYGmyfm0seWYqJEBAMdpK5A0DQ82dYyajskUzUHWItw6GgDIEBcLWRF9o6hM0BZ0uoQcAoA5UCBN5+KcK6+//9TIflwNOjXQc7UkrFgsyvr6uhwdHcnR0VHayQFgKU0FkqnAgeBsOhNrJGEVE9WnAAAYrklEQVQ+zWZTms2m3NzcxNre/avwN169ekVnXwBO0dYPRFOQp42pmp0siCocos6+D8lMIAMAJrjYf8NFpmqZTAUg1Owkh0AGAObgYnOKtnWJTDBVO6QpAKEJcDK3jw4AHGaqYGNdouSx+nVyuOoAwFIUbMnTVMuEyQhkACAFNBMkS9tQcBNBJ0HVZNwRAJACCjbMiyB2MnIFAOAcbc1uBJ3JIZABgBTQ+XM6F5vdNKfNduQsAMxB20ghF9HshnlwJwHAHDTVglBYT6ctUNRUy6QpLSboTBUAOM5EYaK1YFmWpgDNVKGvKQDWlBYT3LwLACAhpgpZ1woTkzQFaJwn/fRcLQBgAU2FrGtNBC7TVMukKS0mrIVhGKadiCRFq2deXl6y+jUANUwEIU8++3KhfVBbEJ+mYFFTWlYhbvlt59EtoFgsyvr6+mh5cABIk62FC9KTlWauZrMpzWZTbm5uYm1PjQwAWCprb+hp0FTrpSktq0CNDAA4joDEHiaCTtf6tpjCXQAAwBSaRqkRuE5GrgDIBJphpiNvpsvCMdqOMwQgE7LSUXIR5E3yaBZKDoEMAAAJo2YnOYxaApAJNJ9MR95kiy3nm1FLAHCLtkJXU2GiLW+QLNeaErl6ASAFrhUmyB4twTiBDAD1tDwwgTSZug80DSk3gbscgHpaHpgi+goTTQg4k2XqPnAtv906GgCquFiwUZhMpyngRPK0BOPu3UkAlmYqADFVsGl5YAJp0nYfaAnGdaQCgCra3qy1PDBF9BUmmpA3ydJ0H2iSmVwpFouyvr4uR0dHcnR0lHZygExwsWCjMJmOvIEJzWZTms2m3NzcxNqeCfEA3ONi3xZtyGNgNibEA7AwCsvkaWu+A2zF0woAMo7aIdiMqxBAJmgrrDX1HzJVO2Qqj7WdK+jGWQeQCdqackwUutoKfFN5rO1cQTcCGQCwFPP0AAQyADKCwno6UzU0pvKYc4V5EMgAUM9EE4qL/SdMFfimmqhM5bGL5wrJsfJq8X1fWq2WPH78WN6+fSvFYlHK5XLayQKQEPpMJIv8hc2sC2R6vZ40Gg3pdDqSy+XE9305ODiQUqkkuVwu7eQBwMoQgABzBDJ/93d/J7lcTp48eZJgch52eHgoL1++HAUtvu/L+fl5qmkCkCz6TCSL/IXNYgcy3W5X1tbW5M/+7M+STM9Mp6en4nmeFAqF0WelUkkuLi5SSxOA5NFnYjJTAQj5C5u9t8yXP/vsM/nxj3889fdffPGFrK+vy+///u/LX//1Xy/zp0REpNVqied5S+8HABb1T//8Lwv9S8I3vvb+Qv8Al8S6op89eybdblc++uijsc/z+bysra1N/d4nn3wilUpF/vIv/1I+++wz+ZM/+ZOlEuv7vpTLZWm32yIiEgSBiIjUarUHv3t1dTX288bGhmxsbCyVHgDZQ78UIBnX19dyfX09+vluuT1NrBqZTz75RP7qr/5KGo2GfPXVV6PPLy8vZTgcTvzOF198IZeXl1Kv10XkXdCzjCho8X1fSqWSVCqVUQBzeHj44Pd3dnZka2tr9O/k5GSp9AAAAHNOTk7GyumdnZ1Y31sLwzCM+0eq1aqcnZ3Jxx9/LJ7nyY9+9CM5ODiQ58+f39t2f39fHj9+LD/5ybu3ly+++EL+9E//NO6fusf3fcnn81Iul6XT6Yw+D4JAHj16JP1+f6zvTCRaBvzNmzdjy4BTIwNgEdqWBQBcMalGZmdnRy4vL8fK77vmurNarZbkcjl59uyZBEEgjUZDCoWCHB8fj9VwvHz5UgaDgQwGg9FnywQxIjLqG1MsFsc+j0Yv9Xq9iYFMZHNzc2ZGAEAcBCRAMhatYJj7jmw0GtJoNO59/r3vfU++973viYjIj370I6lWq/KHf/iHcydoUdOauAAAq0FtFdJg5Op5+vSpPH36VJ49eybD4VCePXu2dA3MJIVCQd6+fTvxd8v2wQEALIeO0EiD0TD4k08+Mbm7e6rV6lj/GJF3fWdEhCUKAADIoLk6+2qQz+el0+mM+sMcHh7K9va2tFqtidtHnX0f6iwEAFgOTUswKW75bd3V0+/3pV6vSy6XkyAIpFgsxppHBgA0cbHQ15w2uMu6Gpl5USMDQKMnn3250PfoT4KscLZGBgDwjou1OsC8uJoBIAUmFnxklBBAIAMAqaBWBDCDOwkALGWiVsckmrqQBq4eALCUqQDAVABCUxfSQCADABlHAAKbEcgAwBxM1V642AyjrakL2aD3jjCsWCzK+vq6HB0dydHRUdrJAWApU7UXmmpBTAUgmoMs2KPZbEqz2ZSbm5tY22fmqnv16hUT4gHABAQg0CSqcIgmxHsIM/sCwBxoWgJWg5l9AeAWU4GDqUBCU0BCUAWbcRUCyARNfVK0MZU3BERIA1cPAMAIgkWkgUAGQCYwNHi6199/mnYSxlCzg3lw1gFkAoXcdPt/8XKh792tSTEVLFKzg3lwZwOApbTVXBAsIg0MvwYASz357MuFvqe9k6629CAdDL8GAMSiLQDQlh7o9l7aCQAAAFgUgQwAALAW9XcAYCmGlAPUyAAAAItlpkamWCzK+vr6aFVNALCdq/OtMGop25rNpjSbTbm5uYm1PcOvAcBSpoZfa+PqcWE+DL8GAMVM1DrQRwYgkAGAVJhoFnK1KYUADfNw8y4AAFjL1QANyeBqAYAUUOsAmEEgAwBzMDWihlqHZDHyKTs4YwAwB1eHPLuG85QdBDIAkHGmai+oBUEamEcGAObgYmFtat4WTfO/uHiesoZ5ZAAgARR0duA8ZQc1MgCQcTQtQSNqZAAAsZgKJAhIkAauOgCAEdTIIA1cPQAAIxjyjDS8l3YCAAAAFpWZGplisSjr6+tydHQkR0dHaScHAJamrSmHZRdgQrPZlGazKTc3N7G2Z9QSAFhK07wtgGlxy2+algAAgLUy07QEAK6hKQcgkAEAazFsGaBpCQAAWIxABgAAWMv6esmDgwPpdrtpJwMArKVtGDcwD6uvwna7Lb1eL+1kAIDVtM3IS2CFeVh71oMgkE6nk3YyAACGaQusoJu1gUy73ZZqtUqNDAAsiWHcsJmVgYzv++J5XtrJAAAnaGuSIbDCPKwctXR2diblcjntZAAAEvCNr72/0D9kk3WBTK/XI4gBAAAiYmHT0mAwkFKpNPf3rq6uxn7e2NiQjY0NU8kCAABLuL6+luvr69HPd8vtaayqkWm321KpVBb67s7OjmxtbY3+nZycGE4dAABY1MnJyVg5vbOzE+t71tTIBEEgIiK5XG6h779582ZsGXBqYwAA0OP4+Fg+/fTT0c9XV1exghlrApleryf9fl+q1eroM9/3RUSkWq1KLpeTRqMx9fubm5tjgQwAwF1MqmefRbt8rIVhGCaQnpWI5pKZdQhXV1eytbUll5eXBDIAkBFPPvtyoe8xqZ4ecctvq/rI3BU1NwEAgGyysg7N931pNBqjWX0PDg7k8PBw4Y7AAAC3mJpUjyYq/axuWoqDpiUAwKJookpPJpqWAABAtlH3BQDAFKz7pB+BDAAAU9DXRT+algAAgLUIZAAAgLUIZAAAgLUIZAAAgLUIZAAAgLUyE8gUi0XZ3d2VZrOZdlIAABnzT//8Lwv9y6Jmsym7u7tSLBZjbc/MvgAAJIwZgucXt/xmgDwAQBXWN8I8OOsAAFV2//wnC33vdu2FtmCIGYKTQyADAHCOiWDIJGqLkkPOAgBUofYC8yCQAQCoYqL2gmAoOwhkAABGaOqXQlNOdnCmAQBGaOuXgmzIzIR4AADAPdTIAACMoF8K0kAgAwAwgn4pSANNSwAAwFoEMgAAwFoEMgAAwFoEMgAAwFqZCWSKxaLs7u5Ks9lMOykAAGCKZrMpu7u7UiwWY22/FoZhmHCaUnV1dSVbW1tyeXkpm5ubaScHAADEELf8zkyNDAAAcA+BDAAAsBaBDAAAsBaBDAAAsBaBDAAAsBaBDAAAsBaBDAAAsBaBDAAAsBaBDAAAsBaBDAAAsBaBDAAAsBaBDAAAsBaBDAAAsFZmAplisSi7u7vSbDbTTgoAAJii2WzK7u6uFIvFWNuvhWEYJpymVMVdBhwAAOgRt/zOTI0MAABwD4EMAACwFoEMAACwFoEMAACwFoEMAACwFoEMAACwFoEMAACw1vtpJ2BevV5Put2uBEEgvu/L4eGhVCqVtJMFAABSYFUgMxgMZDAYSKPREBGRIAjkww8/lH6/L61WK+XUAQCAVbOqaanVakmtVhv9nMvlpNFoSLvdFt/3U0wZAABIg1WBzIsXL+T09HTss/39fRF51+QEAACyxapAZnt7W96+fZt2MgAAgBJW9ZEZDof3PoualKKaGQAAkB1WBTKTtFotKZVKUigUZm53dXU19vPGxoZsbGwkmTQAABDT9fW1XF9fj36+W25PY1XT0l1nZ2fi+750Op0Ht93Z2ZGtra3Rv5OTkxWkEAAAxHFycjJWTu/s7MT63loYhmHCaUtEEASyt7cn3W5XPM+but3V1ZVsbW3JmzdvZHNzc/Q5NTIAAOgxqUZmZ2dHLi8vx8rvu6xtWjo8PHwwiLltc3NzZkYAAID0LFrBYGXTUrValUajMRbEDAaDFFMEAADSYF0g02635fDwcKxzr+/7TIgHAEAGWdW01Ov1pNPpyMHBwVgNTLfbHS1bAAAAssOqzr6PHj2SIAgm/m7aYUSdfR/qLAQAAPSIW35bVSNzcXGRdhIAAIAi1vWRAQAAiBDIAAAAaxHIAAAAaxHIAAAAaxHIAAAAaxHIAAAAaxHIAAAAa2UmkCkWi7K7uyvNZjPtpAAAgCmazabs7u5KsViMtb1VM/sugpl9AQCwT9zyOzM1MgAAwD0EMgAAwFoEMgAAwFoEMgAAwFoEMgAAwFoEMgAAwFoEMgAAwFoEMgAAwFoEMgAAwFoEMgAAwFoEMgAAwFoEMgAAwFoEMgAAwFqZCWSKxaLs7u5Ks9lMOykAAGCKZrMpu7u7UiwWY22/FoZhmHCaUhV3GXAAAKBH3PI7MzUyAADAPQQyAADAWgQyAADAWgQyAADAWgQyAADAWgQyAADAWgQyAADAWgQyAADAWgQyAADAWgQyAADAWgQyAADAWgQyAADAWgQyAADAWgQyAADAWpkJZIrFouzu7kqz2Uw7KQAAYIpmsym7u7tSLBZjbb8WhmGYcJpSdXV1JVtbW3J5eSmbm5tpJwcAAMQQt/zOTI0MAABwD4EMAACwFoEMAACwFoEMAACwFoEMAACwFoEMAACw1vtpJ2Bevu9Lo9GQfD4vuVxOREQqlUrKqQIAAGmwKpAJgkAODg6k3++PgphqtSrtdptgBgCADLKqaenk5ERKpdIoiBERqdfrUq/Xp37n+vp67L8w6/r6Wj7//HPyN0HkcbLI3+SRx8nKev5aNbNvPp+Xer1+r/ZlbW1N+v2+FAqFe9/52c9+Jjs7O/LmzRv57ne/u6qkZgYzJyePPE4W+Zs88jhZruavkzP7+r4vnufd+zyXy0mv10v0b5tao0nTfrStO+Vi3mjKY23HpG0/Jmg6Jk1pMcnFvNGUx1YeU2iJi4uLUETCbrd773ee54W1Wm3i9968eROKSPjmzZul/v4f/MEfLPV9jfsxsY/Ly8tQRMLLy0sV6TG1H01pMZXHmo5J0364hpPfj7Y81pQ3JvajLX9N7SfucVnT2ff8/Hzm74MgmPh5+JuWs5///Odjn29sbMjGxkbsv39zcyNXV1ext7dhPyb2EX1fyzGZ2o+mtJjKY03HpGk/XMPJ70dbHmvKGxP70Za/i+7n+vp6rJ/PV199JSK/LcensaaPjO/7ks/npdvtSqlUGvtdPp+XUqkkrVZr6vcAAIB9Hurjak2NTGRSzcv5+fnYSKbbnjx5IsPhUD744ANZW1sbfT5vjQwAAEjO3RqZMAzl17/+tXznO9+Z+T1rApnt7W0RmdzEFASBPH78eOL33nvvvYkdhAEAgP2sGbWUy+XE87ypfWHuNjcBAAD3WRPIiIiUy2UZDodjnw0GAxGRsTlkDg4O7n3X932pVqtyenoq7XZb2u12sol13KQ83tvbk16vJ0EQSBAE0m635fT0NIXUAQCywprOviLvmpD29vbuLVFwcHAg5XJZRETa7bZUq9WxXs7Tvre3t8fSBguYlMciMtYHSeRdLVm3211l0qy1t7cnjUZD9vf3RUTkxYsXEgSB1Gq10TasM7acOHkcZxtM5/u+tFotefz4sbx9+1aKxeLo2Rz9nmt4OQ/lcSav4aUHeq/YcDgMa7Va2Gq1Rv+NXFxchKVSKbx7WLVaLaxUKvf2k8vlVpJml0zL4zAMw0qlErZarbDRaIT9fj+F1NlLRMb+lUqlsd9fXFyEnueFFxcXo8+i/EY8D+Vx3G0wWbfbDUul0ugaHQ6HY9cs1/DyHsrjMMzmNWxVjcxDTk9PxfM8OTw8HKstWGRpA0w2LY+j3zkd9ScoqiEMgkBKpdK9a7Jer0sQBGNTDPi+L3t7e3JxcbHq5FrpoTyOuw0me/Tokbx8+XKUZ71eTw4PD+Xv//7vJZfLcQ0b8FAei2TzGrZm1NJDpi1fMOt30dIGWTjRJszKYywnn8/PrGI/Ozu7tzhq1Pl9MBhwDcfwUB7H3Qb3RS84t6/DUqk0FqBwDS8nTh6LZPMatqqz7yxnZ2dj7YSRaaOcRN4N6X779m2SyXLKtDyOvH37dtSZ+vT0lI6+BqW5zhjwkFar9eBLDtfwcuLkcVY5USPT6/WmFrCLLm2AcbPyOBJ15ItUq1Wp1+tjn2GyKAjM5XKjazJqpiMYN2NWHs+zDe7zfV/K5fJoNCjXsHkP5XEki9ewEzUyg8GASDVhcfK40+mM/Xx4eCinp6cEizH4vi+1Wk0qlYrUajUZDoejaniCcTNm5fE822BcdP35vi+lUmmUdyLvngEiXMPLipPHkSxew9YHMu12O1Z74LxLG+C34ubxXdHwP6qNH0YQmLw4ecx5mN/tIOX2y06lUpGzs7PRXF9Y3Dx5nMVr2OpAJjoxs4KRRZc2wDtx8ljkXTPStIDF933j6XLdpCCQYNysOIE2wfjDooK1WCyOfR5dl1zDy5snj+/KwjVsdR+ZXq8n/X5fqtXq6LOo0KxWq5LL5aTRaLC0wRLi5nG73Za9vb2x70bBI6MRZqtWq3J4eDjxWoyqkkUIxpfxUB7H3QbzGw6HvFAmLJrxPqvXsNWBTLlcvtcBtd1uS6/XG5urIO7SBrgvbh5HbbK39Xo9yeVyozcCTPZQEMg6Y8uLE2gTjC+uUChM7bAbzeLLNbych/JYJLvXsNVNS5NMulGOj49HawBFWq3WvbZExDMpjw8ODuTs7Gxsm0ajIc+ePaPa+AFxgkCC8eXEyWOC8cVVq9V7fWGiGoDoRYhreDlx8jiz13C6EwubMxwOw0qlEnqeN5qW+fbU17OWNkA8D+Vxt9sdLQdRLpfDbrebYmrt0e12w06nM/o5msp90md3p3e/vQ2mi5PHcbbBdJ7njS1NUi6Xx5aG4Rpe3kN5nNVr2KklCgBb9Xo96Xa7EgSBnJ+fS7VavVfdHi0Wl8/nZTgcZnIGz2XEyeM422CyIAikXq+P5i/J5/P35i/hGl5OnDzO4jVMIAMAAKzlXB8ZAACQHQQyAADAWgQyAADAWgQyAADAWgQyAADAWgQyAJx2cHCQdhIAJIhABoCzTk9PnV4sDwCBDADl6vW67O3tydramuTzeTk8PBxbDmOawWAwmrrd5QXzgKwjkAGgWqPRGK2+3u12pdPp3FvIdBrP8ySXyxHIAA4jkAGgXrfbFc/zxPO8WNsPBoPRQoSe5xHIAA57P+0EAMBDer1e7PVifN+XXq8nr1+/Hn12d9VlAO4gkAGgWhAEEgRB7NFHQRCMLaQ3HA6pkQEcRtMSANWiUUdxamTa7faoSSmSz+cJZACHEcgAUK3b7Uoul5vZP8b3fTk8PJRGozEWtPR6PWm1WjIYDOT09HQVyQWwYmthGIZpJwIApsnn81IoFKTT6aSdFAAKUSMDQK0gCMT3fSkWi2knBYBSBDIA1IpGHk3rH9Nut1eZHAAKEcgAUKvb7YqI3OvAG2FYNQACGQBqzZo/pt1u0+QEgEAGgE5BEIzN0Htbu92Wer0ee6kCAO5iQjwA6lSr1dH8Mb1eb7TWku/78vr1awmCQCqVyr3v1et1abfb8vLlS/E8T87Pz2MvawDATgy/BuAE3/fF933Z39+XFy9eyP7+/tS+NQDcQSADAACsRR8ZAE4KgiDtJABYAQIZAE4YDAbSbrfF930JgkBOTk7SThKAFaBpCYD1ovWVtre35enTpyIi0ul06OgLZACBDAAAsBZNSwAAwFoEMgAAwFoEMgAAwFoEMgAAwFoEMgAAwFoEMgAAwFoEMgAAwFoEMgAAwFoEMgAAwFoEMgAAwFoEMgAAwFoEMgAAwFr/Hx+0jWDoDwBQAAAAAElFTkSuQmCC",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x1358b8550>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = subplots()\n",
    "ax[:plot](ent_spec[:, 1], ent_spec[:, 2], ls=\"none\", marker=\"_\", ms=\"10\", mew=\"1.5\")\n",
    "ax[:set_xlabel](\"\\$L_z^A\\$\")\n",
    "ax[:set_ylabel](\"\\$\\\\xi\\$\")\n",
    "ax[:set_title](\"\\$N=$(N), N_\\\\phi=$(Nphi); N_\\\\mathrm{orb}^A=$(length(subsystemA)), N_e^A=$(NA) : P[0|1]\\$\")\n",
    "ax[:set_xlim](40, 68)\n",
    "ax[:set_ylim](0, 12);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculating entanglement spectrum for NA = 9, LzA = [48.5, 49.5, 50.5, 51.5, 52.5, 53.5, 54.5, 55.5, 56.5, 57.5, 58.5, 59.5, 60.5, 61.5, 62.5, 63.5, 64.5, 65.5, 66.5] ...\n",
      "elapsed time: 686.542615196 seconds\n",
      "\n"
     ]
    }
   ],
   "source": [
    "subsystemA = collect(14:get_Norb(setup))\n",
    "NA = 9\n",
    "LzAvec = collect(48.5:66.5)\n",
    "\n",
    "ent_spec = entanglement_spectrum(psi, setup.states, subsystemA, NA, LzAvec);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjIAAAHUCAYAAAAgOcJbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAPYQAAD2EBqD+naQAAIABJREFUeJzt3c9vI2l62PFHo53VYrOWOOo5ZI0V0FMEDFgIDENqHhIgFze1PuQUmxz/BUMeeJsY5AhYZGaDxcokkCNhk/0XTJN73ENANpAgSAK4m4QPQRAYYHmBXsCXlVTSbIyV13Ll0FMcUeKPEvmS9bxvfT+AMCOKrH7rrbfqffj+3ArDMBQAAAALvZd0AgAAAJZFIAMAAKxFIAMAAKxFIAMAAKxFIAMgceVyWfr9ftLJwAJcJ2hEIAMgUf1+X9rttgRBkHRSMAfXaX36/b588MEH8sEHH8jJycnG//1GozH+98vl8sb//VURyABIVK/XExER3/cTTgnm4Tqt18cffyyXl5fjfN6karUql5eX0ul0Nv5vm0Agg0cbDodSLBZjvy/6qdVqG0idO9rtthSLRTk+Ppbj42NpNBoz3zscDqVcLku5XJZisTj3vZo0Gg0pl8vieZ6MRqOkk7OSefdFEARycnIiw+HQyhYN7dfJ1vL/GIueu77vL2xNifvstg2BzIZEFdLW1pZsbW3NfJgVi0X54IMPZGtrS7LZrKpmviAIpN1uy/Hx8cKHcbvdlk8++UTq9bp0Oh3pdDqSy+WMBTMu5Oc8tVpN9vf3pdPpyGAwkE6nI2dnZ5LNZh+8Nwp46vW6tFot6XQ68vr1a6MPrHXkd/TN3vM88TxPLi4ujKV3k+LcF77vS7/fl+Pj43H+TPtZVRqv0ybKfxzFYlGy2ew4709OTsZf4k5OTuT4+HgczD5G3OduEAQzW8se8+y2UoiNGQwGYT6fD0UkrFarM9/X6XTCQqGwwZQtVigUwnw+H7ZarTCTyYT5fH7meweDQZjJZMLLy8uJ1z3PC4+Ojoylyeb8nGcwGEw9n16vF4pIWCqVxq9dXl6GIhJ2Op2J90av93o9o+kymd93z6NQKBgtG5sS976I/h7dA/d/MplM2Gq1jKQpTddpk+U/jsFg8OAevavVak1Nb6/Xm/qZxz53p/39MceYlQ7tCGQ2qF6vh6PRKMxkMmEmk5n5vlarFQ4Ggw2m7HEW3Qye5029GWa9vixX8vO+Uqn0IAiMZDKZ8O73j1KpFIrI1Pd7njf3Oj2WyfxutVoTlUy1Wp17TBvMuy+q1Wo4Go2m/u3y8tJooJ2m67TJ8h/HrEDlrvv3cBjGCyCWDWQecwxbAxm6ljZoNBqJ53ny8ccfSxAEM6cxjkYjOTo62nDqzOj3++L7/tRm3dFoJK1Wy9i/5Wp+vnnzRj766KOpTcCe54nIN839/X5fMpmMZDKZB+89OjoyOlXWVH4HQSCtVkt6vZ7UajWp1WrWjh15jOja3Ver1aRerxv7d9J0nTZZ/uOIBurm8/mZ79nf3xcRBk2bRCCTgGicyKyHl6YHxWNFo97n3cimuZaf+/v7EgRBrLEI8x6G63pgrprftVpNXr16JfV6ffwTjdGw7VrFdXp6OvX1fr8v2Wx2ZpCzijRcpyTK/zz9fl88z5saWEWi9ETpw+oIZDZkOBzK8fGxiLz7ZhZ9W7j/QPB9f/w+G939BtRoNKTRaIy/zZnkcn5GA3ynVW53B14uEj1M7w4uLJfLUwcML2Iqv9vttpycnDx40N9vaXLNrIqt1WpJtVo19u9wnb4xrfzPs+y9EfF9X4IgmPslrtvtiohIoVCYG+zgcQhkNqTf708U8Ogb2tnZ2dz32cb3fclkMtJut6VUKkm1WpV6vS7ZbFay2ayxb3Iu52cmk5na5B9VSHcrvnldA9Mqm6jr77EV0ar5HXU3lsvlB+tkdLvdcaBbq9VSs3Jsu902PosubdfpseV/nmXvjbufF5G5AWKtVpOjoyNr12tRK+lBOmkxbQCViDwYODdvlsEs0QC8ZX6WGdU/b8BYdE7TZmB4nrfU+U2zzvxcZNP5HTk6Ogo9z5sY2FitVucOdhSRiWsxGo2WGvicZH4vS9N9Mc06Zv9ou07rvgaPLf/zLHtvRAqFQigiUwd1DwaD8WSHaWllsO9qvpVUAAWRUqkk7XZ74tvRMi0WUf+1FkEQyMcff/zg9aOjI2k0GnJ6erqWZtVl87Pb7UqhUIj97ySR341GQ4bDoYxGo4m8Oz09lUajIf1+f+IcfN8Xz/PE9/2JvniTYzFMld910XZf3NVut+XZs2cb+beSvE7rvgaPLf/zrHpvRC0y0yY0PHnyRHq93lrGQoGupY2Y1R99fzCe7/sr9dFqMi1QMdW/bjo/tS3wdZ/v+3J2djb1QZjJZKTX68nZ2dm4cgqCQLrd7vi9qwaNaSy/61av142P3UrjddpE+Y/j7viYu4Ojo59qtUoQs0a0yGzArP7o+4PxhsOhdeM57pv30Hjy5ImIvJtevMp06DTlp8i7FUM7nc7Mc8nn89LpdKTdbo9fq1ar4ynwq37zTyK/o0rJxQGR0QqspmetpO2+iKy7/McRtca4uPy/DWiR2YBZM1BEZDzYr91uy+vXr61a72SaZ8+ezWy2Pj8/F5HVm3BN5mfUBK1VsViU09PThRWP53lSrVbHPyLfnNuqwcCmy2+/35fnz5/L8+fPVz6WRi9fvhQRs918Iul6zty3zvIfR5z1Y7A+BDIbMO9GKpVKIiITTaOPVavVxluwP/bH9MyDaAv6aecSvbbqA9xEfrbbbWm321Kv18X3fanVarG7vDaV37VaTU5OTh6M32m327HSOhwOjUx7X3f5vS9qnl+VpvvirsFgICLmW5s2fZ3iSPIamCr/cUQL82n+UuS0pEcbu240Gi0cNR8ts12v1zeUqtXMG/k+Go1m7nES7Stz36zl+Gcdf9X8zOfz4/RFx7q8vFS1h0yn05mZ/rvL2fd6vTCTyTyYbRHtq3LfY2dmJFV+e72equsRR9xZS9E+SIts+r6w0WPL/zyXl5dLzVqKnnmrbDPBrKXV0CKzZvV6fWGUHjX72tIsOW/VWc/zpFQqPfhGPRwOZTgcyosXLyZej3YDjtsasmp+lstlOTk5efC3TCYzHj+QNN/35ZNPPpEvv/xSjo+PJ36y2exEGqflWxAEUqvVHuS1iIx34Y37rdzF8rsuJlZjjmz6vojS1Wg0pNvtjn+0e2z5nye6xx47GSHKp6g1el3ilq91H0OlpCMpV5VKpfE6BplMZmG0vunNzR6rWq2G+Xx+vOGZiIw3ZZu2QVo+nw9LpVI4GAzCVqsVep43tZWmUCiEnufF+jZpIj/vFvnBYDDxDUxmrEexaUdHR+M8nvZz/9wKhULY6XQm8nrWpnWlUinWGh2m8vvy8jKsVqthp9MJO51OWK1Wx3kctbpEZahUKo2/DUZ/izYtbLVaqtaoiTz2vgjDb9Y3mWfT98VoNJr422g0Upnf0zym/M8T5WVc1Wp14l71PC8sFAozNwedZ1ZLyGPK16wWmcccw9YWGQIZrE30YOl0OgsDhLiLVq0i2hF42r/ZarXUB5PzRF1RrVZrYV4PBoOVFuZ7jKOjo4lgcTQaTVQWnU4n9Dxv3OUVdXv0er0HgWWn07H6GkV6vV7s/N/EfRGG3yxW2Wq1xvesTR5T/jUyEUDE6VraRDqSwPRrrM3R0VGs2RHD4XAjg+Q8zxtvyHh3YORwOJROp2P1suGPWdDv4uJiI/nd7XYlCIKJMhDNIom2sMhkMhPpmfbeSKFQkGKxKMPh0OpZN3G74DZ1X0TTwbUuHhjHY8o/3MMYGSTuyy+/3Nj4imjxrG63K4PBQNrttrx580Z6vZ6Ta5ZMs6kK8vXr11PzdH9/X0aj0fj3x6bFhg0LTdjUfTErPzWt0gzMQyCDRHW7XeMb583jed54YGStVpNSqTSempoGm6yccrnc1Epy1ZVl0zDFdZP3xdHR0Xiwe8T3fRWbQgJxEMggUYVCIZGK6c2bN6moEO97+fLlxO7Z61QoFGR/f3/qLKs4wWO07Huk3W5LPp+3ulsprk3fF51OZ9xS2e12xfd9umtgDcbIACmy6danwWAgZ2dn8ubNGxERGY1G4wXh+v2+1Ov18cJl5XJ5ovJutVrjVoGLiwsZjUbjFVRhVj6fT/30+aS9fPlSXr58Kc+ePVuqnK+yIF+j0ZCzszMRkakb/mq3FYZhmHQiAAAAlqGua2nestL9fn/8ze3k5GRikzAAAJA+qrqWhsOhPH/+fGrzd7QybDRFMAgC+eijj2QwGEir1dp0UgEAgAIqWmR835disShffvnlzK3tW63WxCDFTCYj9Xo99gZ6AADAPSoCGc/zpNPpSL1en7mWx8uXL6XRaEy89uzZMxERpgkCAJBSKgKZOPb39+X8/DzpZAAAAEVUjZGZ5+5KoJGoSylqmQEAAOliTYvMNK1WKzULZAEAgIesaZG5L1p9Mlpca5Z//ud/ll/84hfy/vvvy9bW1vj1nZ0d2dnZWXcyAQBADDc3N3JzczP+PQxD+e1vfytPnz6V996b3e5iZSATBIHUarVYG/394he/WGlfFwAAkJy3b9/KD37wg5l/tzKQKRaL0uv1Yi3H/P7774uIyF//9V/L97///fHrj22RyeVy8vr168cnVvFxTBzj+vpaDg4O5O3bt7K7u5t4ekwdR1NaTOVxlJZ/+Md/Wurz3/32tyaO868+/y9LHed///iPJ46zqlwuJ//tf/yvpT773W9/ayJ//81/Xu446zinVfNY03XSlsem8ybp4/zP//CvjedvEuXmfovMV199JYeHh/I7v/M7cz9nXSBTLpfHuxdHhsPhzHEyUXfS97///bkR3SLb29srV9TajmMqLSIiu7u7atKjKW805XGUlj/47OdLff4Xf/HvJo7z3s53lzpOdA4m8/hffjh9/anH2N3dlf9b//dLfTYK8kyXm1XyODqGlnOK0mUiPf/43reXOsaH3/vO+Bgmy7CW45i83hrKzfX1tYjIxLCQaawKZNrtthSLxYmgxfd98X1/7QN+K5WKc8cxlRZTXMwbTXls+pz+z3/6YyPHWRV5PP8YUcWy6nFMMZGeZz95tdRno2Bc03UycZx/+s0/jP/f1PXWVm7mUbdpZDablXw+/2DbgWin3JOTk4nXe72e1Ov1mYHML3/5y3GT2yotMpju+vpa9vb25Orqyti3Nkwyncerdi1ptMo53c3fP/jpf1/qOFEFaSI9LjJdhp+u2Kpo+jhJX+91PIeTPieR+Oel4q4JgkDOzs7GrSvtdlsuLi4kl8uNtyUoFosSBMHUVXzntcZE42CYobQeOzs78vnnn5O/a2Q6j12sLA//43L9+b/4i3+3ljK8Snru0lCZRFZJi+k8fvOj50aOY4qp670szWV4E9S1yJhGiwFcoqli00TbN2ttLQYmaEqLKdqut4n0aDunVVjVIgMgHk3fkjQFVauOUYiYSpu2FgNNNJUbU0yVPxP3t6lnhKlz2gS9JQOAapq6T0xVcqYq2VUHo0ZMVCbaAgdNwbiptGgOspZl0znZk1IAVn1LisvFik0TF89JG1PBoon728VnxCIEMsAGmHrQafqWlMYHZlymupY0BSGmzklTudHUJSRi5v7W9IzYlPSdMZAATRWSKWl8YMZlqmtJE1PnRLmBaZQoYA5t4wtc5OI4EFNc7GrQdK20DYzVlDc2SffZAwtoe9Bpoqm7TFuLl7ZZVCZo64bRxNR1cjFvNkHPXQI4TFOFZAoP3dk0XW9N4zdcRUtKslKTi7lcTra3t6VSqajamwW6udiSgtm0XW8XK0hNXYnaWpm0lb+kNJtNaTabcnt7G+v9rOwLYCmmKhMTx3l78f+WOsbB/r9Y6nOLaFpd1cVgSMOqs3dpWtnXJazsC2CtNI0L+LeN/7ryMUzS1O2mbbFAzKap3NiEEgYAX9NWWZtYu8XUOWlayZkuGNxFIANgKZoqfW1jHUwxsXaLtnMykR5trTymyp+m8UM2sTflABKlqYLU9hB2scVA0zlpq6xNHdfFpQg2QdfdDwBL0DaLRVNgpW0nbhN5nMbKGrPpudsAWEXTN3TWSpnN1NYCLuaxppl3pqRxlWE9JQqAVUw9sDQFRKbYVAnYSNuYKE0tRJpmE24Kdw1gERMVpLZK1sRxtX0L1VQJmMobU11UJvKYgG82bff3JrAgHmAREwtvaVtMTNOD11TeuLhAmra8MUFb15KJ47hU9lgQD8Baudh6YYqLg4Y10bRhqUma7gVteTOPPSkFYKSC1Da+QBMXAxBtM7pcnLWkKT0ujjlbRM/dBmAhTRWkKZoevNry10QQom22kbY8hv0oUUDKaNupl4ptNk3f9DXRFPxqk8Yyk5onSC6Xk+3tbalUKlKpVJJODrAUDQPw1nnMZWnKF6yftuumqdvNBc1mU5rNptze3sZ6P7OWgA0wVdGamJHgYqWvaSaMSS5Ot9eWHte4lL/MWgIU0dTcq/GBpYWmabTTfjdxzKRpuhdcqvQjmtO2Luk7Y8BiLjY9a6pMXFzp1UVM/cddBDLABjAwdjYTlYmLAZ6IriDPFBenX7vIprKnt7QDDtFcsbjA1X2fTFTYLnaXmaLtemtiU7Cop0QBSCUTlYmLlawpLnaX0cI5m00tKabYm3IAidIUPGiqZEWYjjsPm0aul7Z1ojaB0gBgKaYemJq+QWoKzrQxtfu1pu4ybTSdl/a8usuelAJwkqbBvtpadjTlzbOfvFo5Laa4GESL6Lre2vJmHgIZAEvR1PSs7dujTZXApmkqN9oCVxNMlSGb8sb9uwbAWrg4U0hT94mIrqnKmgbYaiozSB6BDIBEaWqF0NR94ipNg321BUSa0qMpLYvoeYIAgCNcHLPj4jRuTUG0iJn0pHHAuj0pBeAkTRsjauo+McVUdxnsoClQ3BQ9dxuAVDLx4DX18NYUgIiYCUJMdZeZCogIrGCarrt2jXK5nGxvb0ulUpFKpZJ0coDEMKNmNm15o2nMjqm0aDonF9k0tmWWZrMpzWZTbm9vY73f/SfT116/fi27u7tJJwNInLamZxMPXm1jUjQFRC5UbGmhaSB0kqIGh+vra9nb21v4fvvPGIDVTDx4tT28tU15NsFUl5Cmc9JG25cMW+i6+wGsXRpX/oxL2/gNTftQmeoS0nz9XeDifbmIvSkHsJQ0rvwZF+M31i+NFW1cmhZAtIn7JQMAUkzbTtwubhqZxrVbNCHXACzFxbEOLo4DcbFy1NbqoCk9msrepqgr4cPhUL788kup1+sP/ub7vtTrdclms5LJZEREpFQqbTqJAMTNClJb19Kvfv2bpT734fe+Yzglbg5gdpGL9+Uiqs54OBzK8+fPpwYnQRDIycmJDAaDcRBTLpel3W4TzABQRdsAW01c3DRSW3pM0NZ9N4+KQMb3fanVauJ5nuzv7099z9nZmeTz+XEQIyJSq9Xk+PiYQAaAES5WSKZoyhttrQ7a0mOCpu6yRVTkvud50ul0RESk3+9PfU+325Varfbgc0EQyHA4lKOjo7WnE4BO2gZbmqr0NU0Hd7GytqnVAbNZczV83xfP8x68nslkpN/vE8gAKWbTt8fH0LQbsqZK31RaXC03JmhqgVvEikAmCIKZf9vf35fz8/MNpgYA5jNVQWraUFNTpa8pLa6yqdXJipReXFzM/fu8QAeA+2z69ph2mlp2KDdusCKQMeH6+nri952dHdnZ2UkoNUByXOxqMHVMTeck4mZFa6I1xdTYIZtaHdLg5uZGbm5uxr/fr7dnseoqTmt5ubi4mJjJNMvBwcHE759//rl88cUXppIGbISJitbFrgZXx0y4OFXZBBenpePd7OQf//jHj/6cFYFMNCV7WhdTEATy5MmThcd4+/at7O7ujn+nNQY20lbRmqApODPFVIuBthYiE1wMrGDG6empfPrpp+Pfr6+vHzRCTKO3tN+RyWTGU62nyefzC4+xu7s7EcgAaeXi3jummDonUy0Gmgb7akIw5KZlh3xYEciIiBQKBRmNRhOvDYdDERGmXiM1TDzATX1b1/St31TFpumctNEUuHKdcJe60hAEwdSWl9PTUzk+PpYgCMZjYlqt1nghPSANXHyAm+iG0TbYV1OrF0EeXLcVhmGYdCKCIJCzszPxfV+63a6IvGuByeVyUq1Wx+/zfV9arZZks1kZjUaSzWYXbk9wfX0te3t7cnV1RdcSoNDTz36+1OfW0fWhKS2uMhEsujh2CA/Frb9VXNVMJjN1t+v7PM+L9T4AcIGmSl/TNhAujvvB8lQEMgDSS9PATVNp0TQdXNt0e02tKZrSguVxNQAkSlOlYCottBjMZiJvNA08RvL0PEEAYEmufrM2MRBa08BjU7RfN2yWisG+68RgX8B92gbpmgqstJ2XCb/69W+W+tyH3/uO4ZToSgsesmqwLwBooGlAq6tMLBZo6jqx1YEbuNsALEVTd46rYyZMbXXgGm3XCckikAFSRtOMGlNcbQHR1GKgKXA1RdO4HyxPbwkDsBaaAhBtXKzYXAxcTbVUaQ6yEB9XEcBSXKz0tVVsJvJYUwBiiqmWKhdbmdIoNVcjl8vJ9va2VCoVqVQqSScHSAx779iDPF4vF4M8FzSbTWk2m3J7exvr/am5S16/fs30a0DcrBz5Zm0PTRthQqeowSGafr0IdzEA62n7Zq1tfyNMR0DkBko7ABimaX8jUwNjNU1xN5W/BHxu4CoCsB7frGfTNIUbWAcCGQDWc/WbtaYF8TR13xG44i43734AcICJ1hRNwZAprgauWA6lAQAcZqprycVdtOEGAhkAMMzFSt9UK4iJ4zArDHdxNQDga9p2vzZxHE3BkCmaZoUheQQyABKl6VuxixUbrQdwHSUcQKJcDB400RQomqJtbRwkS29JBYANc7FiMxUoahqXYmoAs+ZgDfFxFQEkSlPwQMU2m4vjUlxsrUojrgaARFEprJemQNEUTdslIHmpeYLkcjnZ3t4e76oJANppajHQtLCetuBXU7ebC5rNpjSbTbm9vY31/q0wDMM1pylR0TbgV1dXsru7m3RyAGfw0F2/p5/9fKnP3W0xMHEMk8fRVG5MpUVbHrsibv3NEwXAUmiWh+0Iqt3AVQQApTSNb3FxXIqpFhkXV3K2CYEMgKXw0F0/WgzWy1RQpWkl5zQi1wAshYfu+ploMdDWkuJiAKxp3E8akYsAoJSJ4MHFylJbl5Cm7rI0cq+EAwDUMlHpa+sSQrK4igDgMFOtF5rWkQHuIpABAKVMdH2Yar0wtb+RiXMiqMJdBDIArOfqYEvt6VuGiXMyFVRpG2uD5bh3lwBIHRd3eDaFNU5mY6yNG8h9APiaizs8a1vjxESQp61rSVPgmkbkIgDrudha4CoTQZ6m8Toi+loE0ybdZw8gcZoe3nTDzKbpOpmiLW2aWvJsousqrlEul5Pt7W2pVCpSqVSSTg6Ar7GuyHqZCkCoZGdzMXBNUrPZlGazKbe3t7He795dO8Pr16/nbgMOAKZoqvQ1pcVVpgJgAqJ3ogaH6+tr2dvbW/j+1AQyAHQy8fCmAlg/TQNsXb3eLrYIbsJWGIZh0olYpyiiu7q6okUGwEZoGk9iKi1PP/v5UsdZx4BWTflrkqvntay49bebZw/AGpoe3qbSoqli0ZQWETPpcbW7zNXzWjddJRxA6mh6eGtKC4B4CGQAAAtpGpdiKi2aWgNFdOWxTRgjAyBRmioTTWkRcXM8iab0mBr3g/VgjAwAK2gaw6FpGX4RXWvsmGIiPZqCISTPuqvq+760Wi158uSJnJ+fy5MnT6RarSadLAAY0xY8uMZU/rraRZU2VuViEATSarWkXq+PX+v3+1Iul6XVaiWYMgBJcrUicXGNHW3pMYHANVm67+J7zs7OpFwuT7yWz+elVqsllCIAGmirSExV1iYCLW3Bmon0aNvsEcnSVcIX8H1f+v2+lEqlidf39/cTShEAPKQteHCNtvx1sZXJJrpKwwK5XE7K5bLs7+9LoVAQEZHhcCiZTCbhlAFIEmMdZjN1Ti7mjSlpOEfNrMr9arUqX375pRSLRSkUClIul6XX60mn00k6aQASZKoicbGrwdQ5mTqOiwGRi+dkE+tycTAYyMnJiXS7Xel2uzIYDGJ97vr6euL3nZ0d2dnZWUcSAQAzaJp+zVgbXW5ubuTm5mb8+/16exbrFsSr1WqSzWZFRMYDf1ut1oNxM5FZ24B//vnn8sUXX6wtnQDs4+I3a21dSyYWodO2kJ229Njqiy++kB//+McPXl+0IJ5VgUwUxERBi+/7UiwWZTgcymAwkKOjowefiQKZt2/fTmQELTKADi4GDy7SFBBp2onb5HFM0JSWx5rWInNwcOBWIJPNZmU0Gj14/fj4WPL5/MT6MhG2KAB0c/HbrKbKRFNLiqn0uHhOprh0Pzm3RYHv++J53tS/nZ6eyuvXrzecIiDdND28tdE0ZkJTWkTMpEdbGdKWx2mjqzTM4Xme+L4/9W8XFxeSy+U2nCIg3bQtE4/10nSdtA321cTFc1rEqq6lRqMh5+fnE11Ivu9LrVabOQWbriVgPVxqwjZNU2vVr379m6U+9+H3vmM4Je+YSI+2sufiWBsNnOtaEnm3jky325VyuTxeBO/JkyesIwMkII3f/OLStIv2s5+8WuoY6xoHYio9mrCOUbKsCmRERAqFwnhVXwDJcfVboCaaKjZNaSGIxl08iQDAYdoqfW3p0YS8WY5VY2SWwRgZALbSNGZC0zgQbWNksB5OjpEBgGk0VfgmaUqfi+NAXC03acPVAGA9NjRMF217JFFukkUuAsDXNLUWuOrNj54nnQTjKDfJIpABYD0GSdrDxPRrU4GDi0FVGhHIALCeqSZ6AqLZXOw+MbWmDeUmWXpLmGG5XE62t7elUqlIpVJJOjkAFNJc6SZN05YU2gIHyo1ZzWZTms2m3N7exno/068BwGHadoo2QdNUcKwP068BwHImKlpNLSmmmDonAhI3cBUBQClNs2FcrPRpkXEDVwMADNNUQWpqSRExkzfa1pFBsghkAMAwTd052loPTOSNtnNCsigNAKAUFfZ0plq8WEfGDdwhAYPbAAAYpElEQVQlAGCYtu4c15hq8TK1jgySRSADAIbRkpIumsZEpRHryACwnraKRFt6TNC0doumtIiYW2PHxXKzCtaRAZAa2mafaEuPCZrWbjFVcWsLAFwsN5ug6yoCAMb4hr5eDBp2A11LAKynravhV7/+zVLH+fB735n43cVtATQxlb90La0HXUsAUsPUg5zZMLNpqixdrfC1p08rcg0AlHJxGreL+0e5eJ1sQiADAF8zNdbBVMXm4jd0TQNatQ0adrWlad1Sc/a5XE62t7elUqlIpVJJOjkAFDLVJZT2imXdXG0BMRHkuRAMNZtNaTabcnt7G+v9elK+Zq9fv2awLwAkTNP+US5U+vdpavFaVtTgEA32XUTv1QCADdP2Td/FilZT2kxV+qauk7byZws9JQoAEqapkhVx49v1OmgL8DQtFpjGYEjXXQsAwALaZi1poi0Y34T0nTEArBldDXYwVembuk7aWppske6zB4A7TFUkmroaXKQtwNO2IGPacJcAsJ62AATrRYCHuygNAKynLQDR1mKA6bR15VBulkMgAwBfY0XedNEWAJsoN9qCs02wN+UA8DUCEOAdbWvjbAJ3LYBEmXhgEoBgGXTlzKattWoe7n4AS2GALWznYgCcxuDMvasIYC4CkNlsak4HptG2Ns4mcPcBFjFR0WoLQDQ9MLXlDZAUm4Jze1K6olwuJ9vb2+NdNQEbaapoGWALYB2azaY0m025vb2N9f6tMAzDNacpUdE24FdXV7K7u5t0coCVPP3s50t97m4gQ/fJbOQNoEfc+pu7D7CIiVYQFytdUwGIi3kDuI67FrAIFe10mrrckD605CWLXAQAYAUE0skikAFgPU0zn0TMfUM3cRxNaQHWgRIGwHraKktT39BNHEdTWlylLZBOG113PwAAltEWSKcN068BJEpTl4W2bhi6lpBmTk+/9n1fWq2WPHnyRM7PzyWXy0mhUEg6WUCquLjVgam0mKq8TRxHU1qAdbCuZPb7fanX69LpdCSTyYjv+3JyciL5fF4ymUzSyQNSQ1MAAiC9Ygcyf/M3fyOZTEaePn26xuQsViwW5dWrV+Ogxfd9ubi4SDRNAJanaaCkprQAiCd2INPr9WRra0v+/M//fJ3pmavRaIjneXJ0dDR+LZ/Py+XlZWJpAtLKxb2WNKUFQDwr3bWfffaZfPjhhzODm5/97Gfy8ccfi+d5Uq/X5U/+5E9W+eek1WpNBDEAkkOlv36aBthqSgtwV6wS9uLFC+n1evLxxx9PvJ7NZmVra2vm5z755BMplUryl3/5l/LZZ5+tHMj4vi+FQkHa7baIiARBICIi1Wp14Wevr68nft/Z2ZGdnZ2V0gMA66RpHJKmtMBNNzc3cnNzM/79fr09y3tx3vTJJ5/IX/3VX0m9Xpevvvpq/PrV1ZWMRqOpn/nZz34mV1dXUqvVRORd0LOKKGjxfV/y+byUSqVxAFMsFhd+/uDgQPb29sY/Z2dnK6UHAACYc3Z2NlFPHxwcxPpc7DY/z/Pkj/7oj+Tp06fyZ3/2Z+J5nvz0pz+Vk5OTmQnK5/PjwcH7+/tx/6mp7g7o9Txv/P+lUklqtZoMh8O53U5v376dmIdOawwAxMdAaKzb6empfPrpp+Pfr6+vYwUzj+q8bLVakslk5MWLFxIEgdTrdTk6OpLT09OJFo5Xr17JcDiU4XA4fu1P//RPH/NPPRAFL7lcbuL1aPZSv9+fG8js7u6yIB7gKMZvrB95hXVbdsjHo0tmvV6Xer3+4PUf/vCH8sMf/lBERH76059KuVyWP/zDP3x0gpY1q4sLgPtcHb9BKwiwmJEQ+/nz5/L8+XN58eKFjEYjefHixcotMNMcHR3J+fn51L+tOgYHALQx0QpCaxVcZ7SkfvLJJyYP90C5XJZOpzPxmu/7IiJsUQCkGC0Xs7naWgVEYs1a0qJUKonv+xNjb2q1mpRKpYkBwADS5bvf/tZSPwDsZ92dPBgMpFarSSaTkSAIJJfLxVpHBgDS6M2PniedBGCttsIwDJNOxDrF3QYcAEzRNC7l6Wc/X+pzdC0haXHrb+taZABAO8alAJtDIAMADjM1EFpTKxNwFyUMABxmKpCglQlaEcgAgGFMBwc2h0AGAL5mqvvExe4UgjNo5d7dBgBLovtkNheDM7jBqgXxAAAA7kpNiJ3L5WR7e1sqlYpUKpWkkwNAIbpPgOQ1m01pNptye3sb6/0siAcAANSJW3/TtQQAAKyVmq4lAO4yNduIRd8A+3D3AbCeqdlGLs5aIjiD6yipABJFRbteLgZnwF08CQAkykRFa2q2kbZZSwR5wGKUdgDWM1VxawsANAV5BFXQihIGIFHaWkEwHV1U0IpABkCi+MY+25sfPV/5GAQgcB1PEABQ6tlPXi31uXUEIbScQSsCGQAwTNN4ElMBCC1n0IqSCQCGmerOMdG1pC0A0RTkwQ2UDABQSlPXkimM2YFpBDIAYBjjSYDNSc3u17/3e78n29vbUqlUpFKpJJ0sAA5j76fZXDwnmNVsNqXZbMrt7a387d/+7cLdr1MTyCzKCAAw5elnP1/qc3SfAN+IW3+/t8E0AQAAGEUgAwAArEWnIwAYxmBfYHMIZADAMAamApvD3QYAX2NGDWAf7j4A+JqLi7URnMF1lFQA1qOyns3F4Ay4y/27GIDzTFXWDNIF7EMgAwBfc7GFhuAMrnPvrgWQOlTWs7kYnAF3UcIBWI/KejbGD8F1lFQAcJip8UMERNCKEgYASmkKHpj9BK0IZABAKRPBA+OH4LrUBDK5XE62t7elUqlIpVJJOjkAsBGmWmcIiLApzWZTms2m3N7exnr/VhiG4ZrTlKjr62vZ29uTq6sr2d3dTTo5ABCbpq4lYNPi1t+UdgBQioAEWIy7BAAcRqsOXEdJBQDDNAUP2mYbacobuIGSAQCGaQseNCFvYBqBDAA4jNlGcJ31s5ZOTk6k1+vN/DuzlgBsGt0ns5E3iCsVs5ba7bb0+/2kkwEAE6h0ZyNvYJq1JSoIAul0OkknAwDWhtYLYDFrS3u73ZZyuUyLDABnMTAWWOy9pBOwDN/3xfO8pJMBAAASZmWLTLfblWq1Kt1uN+mkAMDaMOMIWMy6QKbf70uhUEg6GQCwdox1ARaz7i4ZDoeSz+cf/bnr6+uJ33d2dmRnZ8dUsgAAwApubm7k5uZm/Pv9ensWq8bItNttKZVKS3324OBA9vb2xj9nZ2eGUwcA7vqHf/ynpX6AuM7Ozibq6YODg1ifs6ZFJggCERHJZDJLff7t27cTC+rQGgMA8TGDCut2enoqn3766fj36+vrWMGMNYFMv9+XwWAg5XJ5/Jrv+yIiUi6XJZPJSL1en/n53d1dVvYFAECpZYd8WL1FQbSWzLxTYIsCAFgdi/Nh01KxRUHU3QQAWC9NAQlBFe6y8qr6vi/1en28qu/JyYkUi8WlBwIDAOzBeB3cZWUg43metFqtpJMBAAASZmUgAwBIL1Y8xl0EMgAAqzDWBXdZtSAeAADAXQQyAADAWrTPAQA2hqnTMI2SAQDYGKZOwzS6lgAAgLVokQEAbAxTp2FaagKZXC4n29vbUqlUpFKpJJ0cAEglxrpgkWazKc1mU25vb2O93+pNI+Ng00gAWB2DdLFpqdg0EgCwGQzShVYM9gUAANaiRQYAsJCmQbp0c+EurioAYCFTQYCJIIRuLtxFIAMA2BiCEJhGIAMAsIqmbi4kj0AGALAxJoIQTd1cSB5XAwCwMZqCALq53MD0awAAYC09oTEAwDi6T2ZjrI0b3C+pAJBidJ/MloZgLQ3oWgIAANYiHAUAh9F9AtcRyACAw+g+getS07WUy+Xk8PBQms1m0kkBAAAzNJtNOTw8lFwuF+v9W2EYhmtOU6Kur69lb29Prq6uZHd3N+nkAACAGOLW36lpkQEAAO4hkAEAANYikAEAANYikAEAANYikAEAANYikAEAANYikAEAANYikAEAANYikAEAANYikAEAANYikAEAANYikAEAANYikAEAANZKTSCTy+Xk8PBQms1m0kkBAAAzNJtNOTw8lFwuF+v9W2EYhmtOU6LibgMOAAD0iFt/p6ZFBgAAuIdABgAAWItABgAAWItABgAAWItABgAAWItABgAAWItABgAAWOtbSSfgsfr9vvR6PQmCQHzfl2KxKKVSKelkAQCABFgVyAyHQxkOh1Kv10VEJAgC+eijj2QwGEir1Uo4dQAAYNOs6lpqtVpSrVbHv2cyGanX69Jut8X3/QRTBgAAkmBVIPPy5UtpNBoTrz179kxE3nU5AQCAdLEqkNnf35fz8/OkkwEAAJSwaozMaDR68FrUpRS1zAAAgPSwKpCZptVqST6fl6Ojo7nvu76+nvh9Z2dHdnZ21pk0AAAQ083Njdzc3Ix/v19vz2JV19J93W5XfN+XTqez8L0HBweyt7c3/jk7O9tACgEAQBxnZ2cT9fTBwUGsz22FYRiuOW1rEQSBHB8fS6/XE8/zZr7v+vpa9vb25O3bt7K7uzt+nRYZAAD0mNYic3BwIFdXVxP1933Wdi0Vi8WFQcxdu7u7czMCAAAkZ9kGBiu7lsrlstTr9YkgZjgcJpgiAACQBOsCmXa7LcVicWJwr+/7LIgHAEAKWdW11O/3pdPpyMnJyUQLTK/XG29bAAAA0sOqwb4ffPCBBEEw9W+zTiMa7LtosBAAANAjbv1tVYvM5eVl0kkAAACKWDdGBgAAIEIgAwAArEUgAwAArEUgAwAArEUgAwAArEUgAwAArEUgAwAArJWaQCaXy8nh4aE0m82kkwIAAGZoNptyeHgouVwu1vutWtl3GazsCwCAfeLW36lpkQEAAO4hkAEAANYikAEAANYikAEAANYikAEAANYikAEAANYikAEAANYikAEAANYikAEAANYikAEAANYikAEAANYikAEAANYikAEAANZKTSCTy+Xk8PBQms1m0kkBAAAzNJtNOTw8lFwuF+v9W2EYhmtOU6LibgMOAAD0iFt/p6ZFBgAAuIdABgAAWItABgAAWItABgAAWItABgAAWItABgAAWItABgAAWItABgAAWItABgAAWItABgAAWItABgAAWItABgAAWItABgAAWItABgAAWCs1gUwul5PDw0NpNptJJwUAAMzQbDbl8PBQcrlcrPdvhWEYrjlNibq+vpa9vT25urqS3d3dpJMDAABiiFt/p6ZFBgAAuIdABgAAWItABgAAWItABgAAWItABgAAWItABgAAWOtbSSfgsXzfl3q9LtlsVjKZjIiIlEqlhFMFAACSYFUgEwSBnJycyGAwGAcx5XJZ2u02wQwAAClkVdfS2dmZ5PP5cRAjIlKr1aRWq838zM3NzcR/YdbNzY188cUX5O8akcfrRf6uH3m8XmnPX6tW9s1ms1Kr1R60vmxtbclgMJCjo6MHn/nlL38pBwcH8vbtW/nBD36wqaSmBisnrx95vF7k7/qRx+vlav46ubKv7/vied6D1zOZjPT7/bX+26b2aNJ0HG37TrmYN5ryWNs5aTuOCZrOSVNaTHIxbzTlsZXnFFri8vIyFJGw1+s9+JvneWG1Wp36ubdv34YiEr59+3alf//3f//3V/q8xuOYOMbV1VUoIuHV1ZWK9Jg6jqa0mMpjTeek6TiU4fUfR1sea8obE8fRlr+mjhP3vKwZ7HtxcTH370EQTH09/Lrn7O///u8nXt/Z2ZGdnZ3Y//7t7a1cX1/Hfr8NxzFxjOjzWs7J1HE0pcVUHms6J03HoQyv/zja8lhT3pg4jrb8XfY4Nzc3E+N8vvrqKxH5ph6fxZoxMr7vSzablV6vJ/l8fuJv2WxW8vm8tFqtmZ8DAAD2WTTG1ZoWmci0lpeLi4uJmUx3PX36VEajkbz//vuytbU1fv2xLTIAAGB97rfIhGEov/3tb+V3f/d3537OmkBmf39fRKZ3MQVBIE+ePJn6uffee2/qAGEAAGA/a2YtZTIZ8Txv5liY+91NAADAfdYEMiIihUJBRqPRxGvD4VBEZGINmZOTkwef9X1fyuWyNBoNabfb0m6315tYx03L4+PjY+n3+xIEgQRBIO12WxqNRgKpAwCkhTWDfUXedSEdHx8/2KLg5ORECoWCiIi0220pl8sTo5xnfe74+JitDZYwLY9FZGIMksi7VrJer7fJpFnr+PhY6vW6PHv2TEREXr58KUEQSLVaHb+HfcZWEyeP47wHs/m+L61WS548eSLn5+eSy+XGz+bo75Th1SzK41SW4ZUnem/YaDQKq9Vq2Gq1xv+NXF5ehvl8Prx/WtVqNSyVSg+Ok8lkNpJml8zK4zAMw1KpFLZarbBer4eDwSCB1NlLRCZ+8vn8xN8vLy9Dz/PCy8vL8WtRfiOeRXkc9z2Yrtfrhfl8flxGR6PRRJmlDK9uUR6HYTrLsFUtMos0Gg3xPE+KxeJEa8EyWxtgull5HP3N6ah/jaIWwiAIJJ/PPyiTtVpNgiCYWGLA9305Pj6Wy8vLTSfXSovyOO57MN0HH3wgr169GudZv9+XYrEof/d3fyeZTIYybMCiPBZJZxm2ZtbSIrO2L5j3t2hrgzRcaBPm5TFWk81m5zaxd7vdB5ujRoPfh8MhZTiGRXkc9z14KPqCc7cc5vP5iQCFMryaOHksks4ybNVg33m63e5EP2Fk1iwnkXdTus/Pz9eZLKfMyuPI+fn5eDB1o9FgoK9BSe4zBizSarUWfsmhDK8mTh6nlRMtMv1+f2YFu+zWBpg0L48j0UC+SLlcllqtNvEapouCwEwmMy6TUTcdwbgZ8/L4Me/BQ77vS6FQGM8GpQybtyiPI2ksw060yAyHQyLVNYuTx51OZ+L3YrEojUaDYDEG3/elWq1KqVSSarUqo9Fo3AxPMG7GvDx+zHswKSp/vu9LPp8f553Iu2eACGV4VXHyOJLGMmx9INNut2P1Bz52awN8I24e3xdN/6PZeDGCwPWLk8dch8e7G6Tc/bJTKpWk2+2O1/rC8h6Tx2ksw1YHMtGFmReMLLu1Ad6Jk8ci77qRZgUsvu8bT5frpgWBBONmxQm0CcYXiyrWXC438XpULinDq3tMHt+XhjJs9RiZfr8vg8FAyuXy+LWo0iyXy5LJZKRer7O1wQri5nG73Zbj4+OJz0bBI7MR5iuXy1IsFqeWxagpWYRgfBWL8jjue/B4o9GIL5RrFq14n9YybHUgUygUHgxAbbfb0u/3J9YqiLu1AR6Km8dRn+xd/X5fMpnM+BsBplsUBLLP2OriBNoE48s7OjqaOWA3WsWXMryaRXkskt4ybHXX0jTTbpTT09PxHkCRVqv1oC8R8UzL45OTE+l2uxPvqdfr8uLFC5qNF4gTBBKMryZOHhOML69cLj8YCxO1AERfhCjDq4mTx6ktw8kuLGzOaDQKS6VS6HneeFnmu0tfz9vaAPEsyuNerzfeDqJQKIS9Xi/B1Nqj1+uFnU5n/Hu0lPu01+4v7373PZgtTh7HeQ9m8zxvYmuSQqEwsTUMZXh1i/I4rWXYqS0KAFv1+33p9XoSBIFcXFxIuVx+0NwebRaXzWZlNBqlcgXPVcTJ4zjvwXRBEEitVhuvX5LNZh+sX0IZXk2cPE5jGSaQAQAA1nJujAwAAEgPAhkAAGAtAhkAAGAtAhkAAGAtAhkAAGAtAhkATjs5OUk6CQDWiEAGgLMajYbTm+UBIJABoFytVpPj42PZ2tqSbDYrxWJxYjuMWYbD4Xjpdpc3zAPSjkAGgGr1en28+3qv15NOp/NgI9NZPM+TTCZDIAM4jEAGgHq9Xk88zxPP82K9fzgcjjci9DyPQAZw2LeSTgAALNLv92PvF+P7vvT7fXnz5s34tfu7LgNwB4EMANWCIJAgCGLPPgqCYGIjvdFoRIsM4DC6lgCoFs06itMi0263x11KkWw2SyADOIxABoBqvV5PMpnM3PExvu9LsViUer0+EbT0+31ptVoyHA6l0WhsIrkANmwrDMMw6UQAwCzZbFaOjo6k0+kknRQACtEiA0CtIAjE933J5XJJJwWAUgQyANSKZh7NGh/Tbrc3mRwAChHIAFCr1+uJiDwYwBthWjUAAhkAas1bP6bdbtPlBIBABoBOQRBMrNB7V7vdllqtFnurAgDuYkE8AOqUy+Xx+jH9fn+815Lv+/LmzRsJgkBKpdKDz9VqNWm32/Lq1SvxPE8uLi5ib2sAwE5MvwbgBN/3xfd9efbsmbx8+VKePXs2c2wNAHcQyAAAAGsxRgaAk4IgSDoJADaAQAaAE4bDobTbbfF9X4IgkLOzs6STBGAD6FoCYL1of6X9/X15/vy5iIh0Oh0G+gIpQCADAACsRdcSAACwFoEMAACwFoEMAACwFoEMAACwFoEMAACwFoEMAACwFoEMAACwFoEMAACwFoEMAACwFoEMAACwFoEMAACwFoEMAACw1v8HBGVihJAa1EMAAAAASUVORK5CYII=",
      "text/plain": [
       "PyPlot.Figure(PyObject <matplotlib.figure.Figure object at 0x13ec262d0>)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = subplots()\n",
    "ax[:plot](ent_spec[:, 1], ent_spec[:, 2], ls=\"none\", marker=\"_\", ms=\"10\", mew=\"1.5\")\n",
    "ax[:set_xlabel](\"\\$L_z^A\\$\")\n",
    "ax[:set_ylabel](\"\\$\\\\xi\\$\")\n",
    "ax[:set_title](\"\\$N=$(N), N_\\\\phi=$(Nphi); N_\\\\mathrm{orb}^A=$(length(subsystemA)), N_e^A=$(NA) : P[1|1]\\$\")\n",
    "ax[:set_xlim](40, 68)\n",
    "ax[:set_ylim](0, 12);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 0.6.0",
   "language": "julia",
   "name": "julia-0.6"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "0.6.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
